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I.  INTRODUCTION 


There  are  many  problems  in  structural  analysis  which  have  an  axi- 
symmetric  configuration.  An  axisymmetric  problem  has  the  advantage 
that  a  two-dimensional  solution  is  possible.  This  fact  greatly  sim¬ 
plifies  the  analysis  of  such  configurations.  There  are  other  problems 
which  have  an  axisymmetric  appearance,  but  have  nonaxisymmetric  fea¬ 
tures  which  render  a  purely  axisymmetric  solution  inapplicable.  One 
example  of  this  type  of  problem  occurs  in  interior  ballistics  where 
longitudinal  slits  or  slots  are  required  for  sabots.  An  example  of 
this  type  of  configuration  is  given  in  Figure  1. 

One  method  of  analyzing  these  nonaxisymmetric  problems  would  be 
to  use  a  full  three-dimensional  finite  element  approach . (1 *2) *  How¬ 
ever,  if  some  advantage  can  be  taken  of  the  fact  that  these  geometries 
will  have  stress  and  displacement  fields  somewhat  similar  to  the  axi¬ 
symmetric  cases,  some  efficiency  could  be  gained  over  the  complete 
three-dimensional  method  of  analysis.  Consequently,  a  method  that  makes 
use  of  the  two-dimensional  axisymmetric  approach,  but  also  yields  first 
order  nonaxisymmetric  effects,  is  an  appealing  alternative. 

The  objective  of  this  investigation  is  to  develop  a  method  of 
analysis  that  will  take  advantage  of  the  axisymmetric  similarity  but 
still  give  accurate  results  for  stress  and  displacement  fields.  To 
perform  this  analysis  the  structure  is  divided  into  several  segments 
in  the  r-0  plane. (3)  An  axisymmetric  solution  is  obtained  for  each 
segment.  The  segments  are  joined  at  a  certain  number  of  points.  At 
these  points  the  displacements  are  forced  to  match.  Since  the  axisym¬ 
metric  displacements  for  adjacent  segments  won't  match,  a  set  of  per¬ 
turbations  displacements  for  each  segment  is  obtained.  Then  the  axi¬ 
symmetric  and  perturbation  displacements  are  combined  to  give  the 
total  displacements.  Once  the  displacements  are  known,  it  is  a  simple 
matter  to  obtain  the  stresses  and  strains  for  the  body. 

In  addition  to  elastic  loadings,  structural  components  are  often 
loaded  only  once,  but  well  into  the  plastic  region.  Therefore  the 
elastic  analysis  is  expanded  to  include  elastic-plastic  strain  harden¬ 
ing  materials.  For  this  part  of  the  analysis  the  Hill  yield  cri- 
terionH)  for  orthotropic  materials  is  used.  If  the  material  is  iso¬ 
tropic  this  reduces  to  the  von  Mises-Hencky  yield  condition.  To  model 
the  strain  hardening  behavior  the  Prager,  or  kinematic,  hardening 
rule  is  used. 


^Superscripts  refer  to  references. 
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After  the  method  is  developed  a  number  of  numerical  examples  are 
presented.  The  method  is  compared  to  known  solutions.  Also  results 
are  compared  to  the  substructure  analysis  which  resembles  the  per¬ 
turbation  method  presented  in  this  paper. 
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II.  ELASTIC  ANALYSIS 


The  purpose  of  the  present  investigation  is  to  develop  an  approxi¬ 
mate  method  of  analysis  of  elastic  configurations  which  are  almost  axi- 
symmetric  but  have  geometrical  changes  in  the  circumferential  direction 
The  method  of  analysis  which  is  proposed  is  based  on  the  finite  element 
approach . 

The  type  of  geometries  that  this  method  is  designed  to  analyze  is 
represented  by  Figure  1.  If  the  geometry  repeats  at  regular  intervals 
only  one  portion  of  the  total  structure  need  be  analyzed.  Consider 
one  such  portion  as  shown  in  Figure  2.  This  configuration  is  divided 
into  a  number  of  segments  in  the  r-0  plane.  Then  each  segment  is  divi¬ 
ded  into  quadrilateral  elements  as  shown  in  Figure  3.  On  each  face  of 
this  segment  the  finite  element  nodes  are  assumed  to  correspond  to  the 
nodes  in  the  adjacent  segments.  Then  an  axisymmetric  analysis  is  per¬ 
formed  on  each  segment.  The  segments  are  joined  at  four  connecting 
nodes  and  displacements  at  these  nodes  are  forced  to  match. 

Axisymmetric  Analysis 


The  first  step  in  the  analysis  is  to  obtain  an  axisymmetric  solu¬ 
tion  for  each  segment.  The  axisymmetric  solution  is  obtained  from  the 
SAGA  I  Finite  Element  Program  developed  at  the  University  of  Illi¬ 
nois.  This  program  is  capable  of  performing  a  stress  and  displace 

ment  analysis  of  layered,  orthotropic  bodies  of  revolution. 

The  configuration  to  be  analyzed  is  defined  by  its  cross-section 
in  the  r-z  plane.  Then  the  geometry  is  divided  into  quadrilateral 
elements  as  shown  in  Figure  3.  Each  quadrilateral  element  is  further 
divided  into  four  triangular  elements  as  shown  in  Figure  4.  For  each 
triangular  element  the  three  components  of  displacement  are  defined  in 
terms  of  linear  variation  in  the  r  and  z  coordinates  as  follows: 


u  =  a^  +  a^r  +  a^z 

v  =  a4  +  v  +  V 

w  =  ay  +  agr  +  agz 


(2.1) 
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where  u,  v,  w  are  the  displacement  components  in  r,  z,  and  0  di¬ 
rections  respectively.  Using  equation  (2.1)  it  is  possible  to  relate 
the  nine  nodal  displacements  in  the  triangular  element  to  the  unknown 
coefficients  aj^  &2>  etc*  anc*  the  coordinates  of  the  nodes  in  the  r-z 
plane. 

Let  us  now  consider  the  appropriate  strain-displacement  rela¬ 
tions  : 
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By  solving  for  a^  a2 ,  etc.  in  terms  of  the  nodal  displacements ,  it  is 
possible  to  combine  equations  (2.1)  and  (2.2)  and  obtain  the  matrix 
equation: 


{£}  =  [B]  {6} 


(2.3) 


where  {e}  represents  the  six  strain  components,  {6}  is  the  nine  nodal 
displacements,  and  [B]  is  explained  in  Appendix  A. 

Next,  consider  the  stress-strain  relations.  The  axisymmetric 
analysis  allowed  for  an  orthotropic  material.  The  three  axes  of 
orthotropy  are  denoted  by  n,  s,  and  t.  The  basic  stress-strains  re¬ 
lations  are  defined  in  the  n,  s,  t  coordinates  as  follows: 
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(2.4) 


where  En,  Es,  Et,  Gns,  Gst,  Gtn,  vns,  vgt,  and  vtn  are  nine  indepen¬ 
dent  material  coefficients.  The  remaining  three  coefficients  in 
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(E  /E  ) 

sn 

ns 
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(2.5) 


Using  the  angles  a  and  3  as  defined  in  Figure  5  the  stress-strain 
relations  can  be  transformed  into  the  r,  6,  z  system.  The  relations 
given  in  equation  (2.4)  can  also  be  inverted.  This  results  in  the 
matrix  relation: 


{a}  =  [D]  {e} 


(2.6) 
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where  {a}  represents  the  six  stress  components  in  the  r,  0,  z  coordi¬ 
nate  system,  and  [D]  is  further  explained  in  Appendix  A. 

Employing  standard  finite  element  techniques^  the  equilibrium 
equations  for  the  triangular  element  are  obtained  from  the  principle 
of  virtual  work.  These  equations  can  be  written  as: 

[k]  {6}  =  {f }  (2.7) 

where  {&}  is  the  nodal  displacement  matrix  and  { f }  is  nodal  force 
matrix.  In  equation  (2.7)  [k]  is  the  stiffness  matrix  and  is  given  by: 

[k]  =  (  [B]T  [D]  [B]  d  V  (2.8) 

)  V 

where  the  integration  is  over  the  volume  of  the  element. 

By  satisfying  equilibrium  at  each  nodal  point,  a  global  matrix 
equation  can  be  formed  as : 


[K]  {q}  =  {F}  (2.9) 

where  [K]  is  the  global  stiffness  matrix  and  {q}  and  {F}  represent  the 
global  displacement  and  global  force  matrix  respectively. 

Once  this  equation  is  formed  it  can  be  solved  for  the  nodal  dis¬ 
placements  {q}.  Then  the  stresses  and  strains  can  be  found  from 
equations  (2.3)  and  (2.6). 


Nonaxisymmetric  Analysis 


While  the  axisymmetric  solution  for  each  segment  is  being  deter¬ 
mined,  the  nonaxisymmetric  force  and  stiffness  matrices  are  calculated 
and  stored.  In  assembling  the  nonaxisymmetric  stiffness  matrix  the 
following  displacement  functions  are  used  for  each  triangular  element: 


u 
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al  +  V  +  a3Z  + 


a.  +  a_r  +  a^  z  + 
4  5  6 


a7  +  V  +  V  + 
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0(b4+b5r+b6z) 
0(b7+bgr+bgz) 


(2.10) 
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Also  the  nonaxisymmetric  strain  displacement  equations  are  used. 
These  are: 
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(2.11) 


where  again  u,  v  and  w  are  the  displacements  in  the  r,  z  and  0 
directions . 

These  equations  along  with  equation  (2.6)  can  be  combined  to  form 
a  nonaxisymmetric  stiffness  matrix  as  follows: 

[S1]  =  [bi]T  [d]  [bi]  d  v  C2-12) 

where  [SjJ  is  the  nonaxisymmetric  stiffness  matrix  and  [B^]  is  the 
nonaxisymmetric  strain-displacement  matrix  obtained  by  combining 
equations  (2.10)  and  (2.11).  The  matrix  [BjJ  is  defined  as: 


{£}  =  [BJ  {6}  C2-13) 

where  {6}  is  the  eighteen  nonaxisymmetric  displacement  components 
(i.e.,  3  displacements  u,  v,  and  w  at  each  node  on  both  faces  of  the 
segment.)  The  [BjJ  matrix  is  derived  in  Appendix  A. 
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In  general,  the  displacements  generated  by  the  axisymmetric 
solution  will  not  be  compatible  with  the  displacements  generated  for 
an  adjacent  segment.  Consequently,  some  changes  need  to  be  made  in 
these  displacements  to  obtain  the  real  situation.  The  total  displace¬ 
ment  is  now  represented  by  two  components.  If  we  call  the  total  dis¬ 
placement  {u},  we  can  write: 


{u}  =  {u  }  +  {u  } 

a  p 


(2.14) 


where  {ua}  is  the  displacement  from  the  axisymmetric  solution  and 
{up}  is  the  additional  perturbation  displacement. 

For  each  side  of  a  segment  the  perturbation  displacements  are 
expressed  as  follows: 


u 

v 

w 


a2r  +  Z(a3 
a^r  +  z(a7 

a10r  +  z(a 


+ 
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a4r) 


*  a12r) 


(2.15) 


where  r,  z,  and  0  refer  again  to  coordinates  in  Figure  1.  Similar 
expressions  can  be  written  for  the  other  side  of  the  segment. 

If  the  total  displacements  are  forced  to  match  for  two  adjacent 
segments  at  four  nodal  points,  a  sufficient  number  of  equations  to  be 
able  to  solve  for  the  coefficients  of  the  perturbation  displacement 
equations  will  be  obtained.  These  nodes,  where  the  solutions  match 
are  called  connecting  nodes.  Then  for  each  segment: 


{u  }  =  [R]  {p}  (2.16) 

where  {ucp}  is  the  perturbation  displacement  at  the  connecting  nodes, 
[R]  is  a  matrix  of  the  connecting  nodal  coordinates  further  defined 
in  Appendix  A,  and  {p }  is  the  unknown  perturbation  coefficients. 
Solving  for  {p}  yields: 

{p}  =  [R]"1  {ucp}  (2.17) 
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Also,  for  each  quadrilateral  element  in  a  particular  segment: 


=  [G]  {p}  (2.18) 

where  [G]  is  the  matrix  of  nodal  coordinates  for  the  element  further 
defined  in  Appendix  A.  Substituting  equation  (2.17)  into  equation 
(2.18)  yields: 

{up>  =  [G]  [R]"1  {ucp}  (2.19) 

Also  for  each  element: 

=  [sp  ({ua>  +  {Up})  (2.20) 

where  [SjJ  is  the  nonaxisymmetric  stiffness  matrix  and  {f }  is  the  in¬ 
ternal  force  vector.  Now,  applying  the  principle  of  virtual  work  to 
equation  (2.20)  gives: 


6uJ  {f}  =  <5 u^  [Sp  ({u&}  +  {Up}) 

From  equation  (2.19),  it  can  be  seen  that: 


(2.21) 


(2.22) 


Substituting  this  into  equation  (2.21)  and  summing  the  contribution  of 
all  the  elements  in  the  segment,  results  in: 


WI  =  <5^p2([R]"1T[G]T[S1]{ua} 

+  [R]"1T[G]T[S1][G][R]-1{ucp}) 


(2.23) 


Also,  there  are  contributions  from  the  external  body  forces.  This  can 
be  written  as: 
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(2.24) 


6uT  {f  }  =  6uT  {P  } 

p  E  pi 

where  { P ^ }  is  the  external  body  force  vector  for  an  individual  element. 
The  same  process  of  substitution  and  summing  over  all  the  elements 
gives : 


W  =  6uT  E  [R] -1T[G]T{P  }  (2.25) 

E  cp  L  J  L  J  1 

Also  there  will  be  a  contribution  from  the  concentrated  loads.  This 
contribution  can  be  stated  as: 

6u*  {fc}  =  6u^  {g}  (2.26) 

p  p 

where  {g }  is  the  concentrated  force  vector  for  the  element.  Then,  for 
the  whole  segment: 

Wc  =  6uJp  E  [R]"lT[A]T{g}  (2.27) 

where  [A]  is  a  matrix  relating  the  perturbation  displacements  at  the 
nodes  where  the  concentrated  forces  are  applied  to  the  perturbation 
coefficient  matrix.  That  is: 


{up}  =  [A]  {p}  (2*28) 

at  nodes  where  concentrated  forces  are  applied.  The  [A]  matrix  is 
presented  in  more  detail  in  Appendix  A. 

Balancing  the  virtual  work  done  by  all  the  forces  gives: 


W-  +  W  =  WT  (2.29) 

E  c  I 


Here,  the  following  substitutions  can  be  made: 
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(2.30) 


where  the  [R]  has  been  brought  outside  the  summation  since  it  is 
constant  for  the  entire  segment,  and  the  summation  is  over  all  the 
elements  in  the  segment.  The  terms  in  the  summation  can  be  easily 
calculated  during  the  axisymmetric  solution  and  later  be  multiplied 
by  the  [R]“l  terms.  Then  {F^},  [SjJ  ,  and  t F j }  can  be  stored  for  each 
segment.  These  substitutions  give: 


[s  ]  {u  }  + 
k  cp 


{v  - 


{V 


(2.31) 


Since  it  is  the  total  displacements  that  must  match  at  the  connecting 
nodes,  equation  (2.31)  can  be  written  as: 


[sk]  <ucT)  -  [Sk]  {uca}  *  {Fj)  =  (Fe) 


(2.32) 


where  {uca}  is  the  axisymmetric  displacement  at  the  connecting  nodes 
and  {ucTj:  is  the  total  displacement  at  the  connecting  nodes.  Every¬ 
thing  is  known  except  {ucT).  Combining  terms  in  equation  (2.32) 
yields: 


[sk]  (ucT)  =  {Fc}  (2.33) 

where  {Fq}  is  a  combination  of  {F^},  {F j } ,  and  [Sj(]  {uca}.  There  will 
be  one  matrix  equation  for  each  segment.  Now  [S^\  can  be  assembled 
into  a  banded  global  stiffness  matrix  for  the  overall  structure.  Then 
appropriate  boundary  conditions  can  be  applied  to  restrict  the  rigid 
body  motion,  preserve  symmetries,  or  comply  with  external  restraints. 
Then  equation  (2.33)  can  be  solved  for  the  {ucj}  vector,  giving  the 
total  displacements  at  all  of  the  connecting  nodes. 
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Now  it  is  a  simple  matter  to  work  back  from  this  point  to  obtain 
the  stresses  and  strains  for  each  element  in  each  segment.  For  a 
given  segment  {ucp}  and  {uca)  are  known,  and  therefore,  {ucp}  can  be 
calculated  from: 


} 


{ucT}  '  fuca) 


(2- 34) 


Then  knowing  {ucp}  and  [R]  ,  {p}  can  be  obtained  from  equation 

(2.17).  Then  for  each  element  the  perturbation  displacements  {u  } 
can  be  obtained  from  equation  (2.18).  Having  determined  (ua)  an! 
(up),  it  is  easy  to  get  the  total  displacement  from  equation  (2.14). 


At  this  point,  it  is  simple  to  apply  the  basic  finite  element 
relations  developed  earlier: 


{e}  = 

[Bx] 

{6} 

{a}  = 

[D] 

{e} 

(2.35) 


Comparison  to  Substructure  Analysis 

This  method  of  analyzing  nonaxisymmetric  structures  is  geometric¬ 
ally  similar  to  the  substructure  analysis.  The  basic  difference  be¬ 
tween  the  substructure  analysis  and  the  perturbation  analysis  is  in 
the  way  total  displacements  are  calculated.  In  the  perturbation 
analysis  the  displacements  are  the  sum  of  the  axisymmetric  and  the 
perturbation  displacements.  In  an  equivalent  substructure  analysis 
the  internal  degrees  of  freedom  are  eliminated  in  each  substructure 
and  a  set  of  equations  similar  to  equation  (2.15)  is  obtained  except 
instead  of  representing  only  part  of  the  displacements,  they  represent 
the  total  displacements. 

With  these  differences  the  perturbation  analysis  outlined  pre¬ 
viously  can  be  converted  into  a  substructure  analysis  by  modifying 
equation  (2.19)  to  read: 


(u)  =  [G]  [R]'1  {ucT> 


(2.36) 
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where  {u}  is  the  displacement  matrix  for  an  element  and  [G]  and  [R] 
are  defined  as  before.  It  is  possible  to  solve  for  {uct)  by  using 
equation  (2.33)  by  replacing  {Fc}  with  {Fp}.  That  is: 


[Sk]  {ucT)  =  {Fe} 


(2.37) 


with  {Fp}  and  [SjJ  defined  as  in  equation  (2.30).  Although  there  is 
major  conceptual  differences  between  the  perturbation  and  substructure 
methods ,  there  are  only  minor  changes  required  in  the  computer  pro¬ 
grams.  Therefore *  it  is  possible  to  obtain  solutions  for  both  methods 
and  compare  answers . 
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III.  NUMERICAL  RESULTS  OF  ELASTIC  ANALYSIS 


To  test  the  accuracy  of  this  analysis  three  examples  were  used. 

On  each  example  both  a  perturbation  and  a  substructure  analysis  were 
performed.  Whenever  possible  the  results  were  compared  to  known 
solutions . 

Axisymmetric  Problem 

The  first  test  of  the  perturbation  method  is  that  if  the  confi¬ 
guration  is  indeed  completely  axisymmetric,  the  solution  should  be 
identical  to  the  axisymmetric  solution.  To  test  this  requirement,  an 
example  of  an  axisymmetric  disk  with  a  center  hole,  loaded  with  in¬ 
ternal  pressure,  was  used.  Both  the  perturbation  and  the  substructure 
analyses  were  performed. 

In  both  analyses  two  segments  in  the  circumferential  direction 
were  used.  Each  segment  was  divided  into  four  elements  in  the  radial 
direction  and  one  element  in  the  axial  direction.  The  disk  had  an 
inner  radius  of  5  inches  and  an  outer  radius  of  15  inches,  and  was  2 
inches  thick.  Thus,  each  element  measured  2.5  inches  by  2  inches. 

In  Table  1  the  results  for  the  radial  stress,  arr,  and  the  cir¬ 
cumferential  stress,  cfgg,  are  presented.  The  axisymmetric  results 
presented  were  obtained  from  the  SAGA  I  finite  element  program  (6,7), 
These  results  agree  well  with  analytic  results.  From  this  table  it 
can  be  seen  that  the  perturbation  analysis  does  yield  answers  that 
agree  well  with  the  axisymmetric  solution  for  both  arr  and  Oqq.  It 
can  also  be  seen  that  the  substructure  analysis  gives  reasonable  ans¬ 
wers  for  egg  but  shows  large  errors  in  computing  the  radial  stress, 

arr- 

Nonaxisymmetric  Disk 


The  second  example  tested  was  a  nonaxisymmetric  disk  under  inter¬ 
nal  pressure  as  shown  in  Figure  6.  Again  both  a  perturbation  and  a 
substructure  analysis  were  performed.  For  this  configuration  six 
segments  in  the  circumferential  direction  were  used.  Three  of  these 
segments  were  the  same  as  those  described  in  the  previous  section, 
and  the  other  three  were  of  the  same  thickness  but  had  only  two  ele¬ 
ments  in  the  radial  direction.  These  elements  had  an  outer  radius  of 
ten  inches. 

Since  this  problem  is  a  plane  stress  problem  in  the  r-0  plane, 
the  solution  obtained  by  the  perturbation  analysis  can  be  compared 
to  a  plane  stress  solution.  A  plane  stress  solution  was  obtained 
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using  the  SAAS  Finite  Element  Program.  In  Figure  7  results  for 
radial  displacement  as  a  function  of  radial  displacement  are  presented. 
Results  for  circumferential  stress  and  radial  stress  as  a  function  of 
radial  location  are  presented  in  Figures  8  and  9  and  Tables  2  and  3. 

It  can  be  seen  that  the  perturbation  solution  agrees  well  with  the 
plane  stress  solution  which  is  known  to  be  accurate.  Furthermore, 
the  substructure  results  disagree  with  the  plane  stress  results  for 
radial  stresses  and  radial  displacements.  The  circumferential  stresses 
do  agree  well  for  the  substructure  analysis.  These  results  are  not 
shown  in  Figure  8  for  the  sake  of  clarity  since  the  results  of  all 
three  analyses  are  very  close. 

Three  Dimensional  Tube 


The  third  example  considered  is  a  three  dimensional  tube  shown 
schematically  in  Figure  10.  Again  six  segments  in  the  0-direction 
were  used.  Four  elements  were  used  in  the  axial  or  z-direction.  The 
end  labelled  the  "unsymmetric  end”  in  Figure  10  has  an  end  view  the 
same  as  shown  in  Figure  6.  The  ,faxisymmetric  end”  has  an  end  view  as 
shown  in  Figure  11.  The  tube  measured  eight  inches  in  the  axial 
direction.  The  geometric  discontinuity  is  at  the  mid-point  of  the 
tube.  As  before,  the  tube  is  subjected  to  an  internal  pressure. 

Once  again  the  problem  was  solved  using  both  the  substructure  and 
the  perturbation  methods.  Results  for  the  circumferential  stresses 
are  presented  in  Table  4,  and  for  the  radial  stresses  in  Table  5. 

Since  this  problem  is  a  true  three  dimensional  problem,  no  analytic 
solution  was  available.  However,  each  end  should  approach  a  limit 
solution  as  represented  by  the  previous  two  examples.  From  Table  4 
and  5  it  can  be  seen  that  the  stresses  do  approach  the  so-called  limit 
solutions.  Also  it  can  be  seen  that  again  the  substructure  method  does 
not  predict  the  radial  stresses  well.  A  boundary  condition  of  the 
problem  is  that  the  outer  boundary  is  stress  free.  The  perturbation 
method  fulfills  this  boundary  condition  while  from  Table  5  it  is 
obvious  that  the  substructure  method  violates  it. 

Effect  of  the  Choice  of  Connecting  Nodes 

One  ambiguity  in  this  method  is  the  choice  of  the  connecting 
nodes.  If  the  accuracy  of  the  solution  depends  on  this  choice,  the 
same  constraints  must  be  placed  on  this  choice.  To  see  what  effect 
the  choice  of  connecting  nodes  has  example  three,  the  nonaxisymmetric 
tube,  was  solved  using  three  sets  of  connecting  nodes.  These  three 
sets  of  nodes  are  shown  in  Figure  12,  where  the  circles  represent  the 
connecting  nodes.  The  choices  ranged  from  widely  spaced  as  shown  in 
Figure  12a  to  very  close  together,  in  fact  all  nodes  of  one  element. 
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as  shown  in  Figure  12c.  In  Table  6  results  are  shown  for  these  three 
cases.  Results  for  the  other  stresses  show  similar  changes.  It  can 
be  seen  that  the  results  are  insensitive  to  the  choice  of  connecting 
nodes.  Thus,  this  choice  can  be  arbitrary  and  no  restraints  need  be 
placed  on  the  generality  of  the  method. 
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iv.  the  elastic-plastic:  analysis 


Many  applications  of  the  perturbation  analysis  described  in 
Chapter  II  occur  in  the  field  of  interior  ballistics.  In  many  cases 
the  structures  involved  must  survive  an  extreme  environment  only 
once.  Therefore,  the  components  are  often  loaded  well  into  the  plas¬ 
tic  range.  To  analyze  such  problems  the  perturbation  method  was  ex¬ 
panded  to  include  elastic-plastic  strain  hardening  material  behavior. 

Incremental  Equation  for  a  General  Strain  Hardening  Material 

To  describe  the  plastic  work  hardening  behavior  of  a  material 
three  things  are  needed: 

a)  an  initial  yield  condition 

b)  a  flow  rule  relating  plastic  strain  increments 
to  the  stress  and  the  stress  increment 

c)  a  hardening  rule 

f41 

For  this  work  the  Hill's  yield  criterion^  J  for  orthotropic 
materials  was  chosen.  This  reduces  to  the  von  Mises-Hencky  yield 
condition  if  the  material  is  isotropic.  The  Hill's  condition  will  be 
discussed  more  fully  later,  but  it  can  be  represented  symbolically  as: 


F 


(a.  .) 

ij 


=  1 


(4.1) 


at  the  point  where  yielding  occurs. 

The  Prager  or  the  kinematic  strain  hardening  rule  was  used^. 
This  hardening  rule  assumes  that  the  yield  surface  retains  its  initial 
size  and  shape  but  undergoes  a  translation  in  the  direction  of  the 
plastic  strain  increment.  This  is  illustrated  in  Figure  13.  After 
plastic  flow  begins  the  yield  surface  can  be  written  as: 


F  (a.  .  -  a.  .)  =  1  (4-2) 

ij  i  Y 

where  the  ot-jj's  represent  the  translation  of  yield  surface.  The 
i  assumption  that  this  translation  is  in  the  direction  of  the  plastic 

strain  increment  can  be  written  as: 
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da. .  =  c  de?.  (4.3) 

ij  ij 

where  c  is  a  constant  for  the  material. 

The  flow  rule  chosen  was  also  due  to  von  Mises,  namely: 


D 

d eii  =  3^  d  *  ;  d  X  >  0  (4.4) 

ij 

The  scalar  dX  in  equation  (4.4)  can  be  determined  by  the  fact  that  the 
stresses  remain  on  the  yield  surface  during  plastic  flow.  This  condi¬ 
tion  can  be  expressed  as: 


(da.  .  -  da.  .)  =  q 

ij  lv  9a.  . 

1 3 


(4.5) 


Substituting  from  equations  (4. 3) -(4. 4)  into  equation  (4.5)  gives: 


9F  A  9F 

c  9a. .  d  A)  9a. . 
iJ  ij 


Solving  equation  (4.6)  for  dX  gives: 


(4.6) 


dX  .  [if  da..}  t|f  f  }-1 
c  3a..  ij  3akl  3akl 


(4.7) 


The  total  strain  increment  is  the  sum  of  the  elastic  and  the  plastic 
strain  increments.  Thus: 


j  e 
de.  . 
ij 


de. .  -  de?. 


ij 


ij 


(4.8) 


where  de^_.  is  the  total  strain  increment. 


gives : 


0 

Substituting  for  de^.  in  the  elastic  stress-strain  relationship 


da.  . 

ij 


Eijkl  ^d£kl  "  dekd 


(4.9) 


where  E. is  the  elastic-stress-strain  tensor. 
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Equation  (4.5)  can  be  rewritten  as: 


(da.  .  -  c  de? . )  =0 

i  J  13  ^aij 


(4.10) 


Then  substituting  equations  (4.4),  (4.7)  and  (4.9)  into  equation 
(4.10)  gives: 


{E.  .  [de  -  fj-  dX] 

ljkl  kl  9a^ 


9F  a\\  n 

"  c  to.  .  dX}  to.  .  =  ° 

13  13 


(4.11) 


This  can  be  rewritten  as: 


E.  de  |J 
ijkl  kl  9a. . 

J  i3 

9F  9F  ,  .. 

*  c  to.,  a?..1  dX 
13  13 


(F  9F  9F 

iikl  9a.  .  9a,  . 
J  13  kl 


(4.12) 


Then  defining  the  bracketed  term  as  1/D  where  D  is  the  scalar: 

„  9F  9F  9F  9F  t  -1 

D  "  {Eijkl  9a,  .  9a..  +  c  9a..  9a.!  (4.13) 

J  kl  13  13  13 

equation  (4.12)  can  be  written  as: 


dX  =  D  E.  ...  ■?—  de.  . 

13 kl  9a. .  kl 

13 


(4.14) 


Substituting  equation  (4.14)  into  equation  (4.4)  gives  a  relation  be¬ 
tween  the  total  strain  increment  and  the  plastic  strain  increment  as: 


de? . 

13 


9F  9F 

mnkl  9a  9a. .  dekl 
mn  13 


(4.15) 
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or : 


C.  ...  d 
ljkl 


'kl 


where  is  defined  from  equation  (4.15). 

Then  equation  (4.9)  can  be  written  as: 

da  =  [E.  ..  -  E.  ,  C  .  .1  d  e.  . 

kl  ljkl  klmn  mni3J  13 

or : 


da.  .  =  A.  .  d  e.  . 

13  13k!  kl 


where : 


1  =  TF  -  F  C  1 

ijkl  ijkl  ijmn  mnklJ 


(4.16) 


(4.17) 


(4.18) 


(4.19) 


For  the  isotropic  case,  c  can  be  evaluated  from  a  uniaxial  test. 
It  can  be  shown  that : 


2  T 

C  3  E-E?  (4.20) 

where  ET  is  the  tangent  modulus  as  shown  in  Figure  14,  and  E  is  Young's 
Modulus.  In  general  ET  will  depend  on  the  stress  state  at  a  given 
time.  Thus  by  knowing  the  elastic  constants  and  the  tangent  modulus 
the  stresses  can  be  obtained  from  the  total  strains. 

Application  of  General  Strain  Hardening  Relations  to  an  Orthotropic 

Material  —  —  c 


In  applying  the  equations  developed  in  the  previous  section  to  an 
orthotropic  material,  the  Hill's  yield  condition  was  used.  This  yield 
condition  can  be  expressed  as: 


24 


F(a.  .) 
ij 


11 


ll 


a 


22 


a 


33 


22 


33 


12 


12 


°23 

+  -44-  + 


'23 


13 


13 


+  Yll°22a33 


+  Y22aiia33  +  Y33aHa22  1 


(4.21) 


where  ,  Y22  and  Y33  are  the  yield  stresses  in  simple  tension  in  the 
three  orthotropic  directions  and  Yj^,  ^23_an<^  ^13  are  correspond¬ 
ing  yield  stresses  in  simple  shear.  The  Y  terms  are  functions  of  the 
yield  stresses.  These  functions  are: 


Y 


11 


Y 


22 


(4.22) 


After  strain  hardening  has  occurred  equation  (4.21)  becomes: 

.2 


(airaiP  ♦ 


(^a22"a22')  +  ^33  a33^ 


11 


22 


'33 


(-ai2  ai2')  +  ^a23"°23^  +  ('ai3"ai3') 


12 


23 


13 


Yll(-a22_a22')  ^a33  a33->  +  Y22  ^33  ^ll"01!!-1 


(4.23) 


+  Y33(airail)(a22-a22) 


1 
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Then  evaluating  the  terms  gives: 

ij 


IL  .  2(011-a11) 
8  a 


11  y2 
11 


— — —  +  Y22(-a33  +  Y33(-a22_a22'1 


3  F 
3a 


22 


22 — —  +  Y33('ail"ail'1  +  Yll('a33“a33') 
Y22 


3_F 
3  a 


33 


— — ^ — —  +  Yll(-a22  °t22')  +  Y22(-ail"ail'1 
Y33 


3  F 
3a 


12 


32 

3  a 


23 


2C°12-ai2) 


12 


2  (-a23_a23') 


23 


32 

3  a 


13 


13 


(4.24) 


There  still  remains  the  question  of  evaluating  the  hardening  para¬ 
meter,  c.  This  was  evaluated  from  a  uniaxial  test  in  the  isotropic 
case,  but  there  are  three  uniaxial  tests  which  will  give  three  dif¬ 
ferent  results  in  the  orthotropic  case.  To  solve  for  the  hardening 
parameter  two  assumptions  were  used  which  reduce  the  generality  of 
the  analysis  but  increase  its  simplicity.  These  assumptions  are: 

a)  transverse  isotropy  (i.e.,  Y22  =  Y33) 

b)  there  is  one  predominant  direction 

Then  by  taking  a  uniaxial  tension  loading  in  the  preferred  direction, 
such  that  a  ^  0  and  all  other  stresses  are  zero,  equation  (4.24) 
becomes : 
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a_F 
3  a 


11 


3F 

3a 


22 


2^ail~aiP  -  a  Y  -  a  Y 
- ^ — —  33  22  22  33 


11 

2a 


—  "  a33  Y11  "  Y33  (airaiP 


22 


3F 
3  a 


33 


3_F 

3a 


12 


2a 


—  "  “22  Y11  +  Y22  ^°11 "“lP 


33 


3^ 

3a 


13 


3_F 

3a 


23 


(4.25) 


By  assuming  transverse  isotropy  equation  (4.22)  becomes: 
a-  1  2 


11 


11 


22 


22 


33 


11 

1 

fll 


(4.26) 


Also  as  a  result  of  transverse  isotropy  in  this  uniaxial  loading  case: 


a22  “33 


(4.27) 


Then,  since  the  plastic  flow  is  incompressible  and  da. .  is  proportion¬ 
al  to  de?.,  it  follows  that:  ^ 

i} 


a  +  a  +  a 
11  22  33 


0 


(4.28) 


and  substituting  from  equation  (4.27)  gives: 


11 


-2a 


22 


-2a 


33 


(4.29) 
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By  substituting  equations  (4.26)  and  (4.29)  into  equation  (4.25): 


9F 

90 


11 


9_F 
9  0 


22 


9F 

90 


33 


2°11  3all 


11 

1  2°11  -  3ail 

2  v2 


11 

20,  ,  -  3a 


1  11 
2 


11 


11 


Noting  the  common  factor  (20^  -  3a^p  gives: 


(4.30) 


1 

2  90 


11 


9F 

90 


22 


9F 

90 


33 


Evaluating  dA  from  equation  (4.7)  gives: 
1/c  do, 


dA 


11 


3/2  9F/90 


11 


(4.31) 


(4.32) 


Then  from  equation  (4.4): 


or: 


I  I 

3  c 


do 


11 


(4.33) 


doll 

= 

3  C 

< 

2 

This  result  is 
now  becomes: 

identical  with  the  isotropic  result. 

r 

2 

En‘EnT 

Li  — 

3 

EirEnT 

(4.34) 

Equation  (4.21) 


(4.35) 
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t 

By  assuming  that  this  direction  is  the  predominant  one  the  value  for 
"c"  given  by  equation  (4.35)  can  be  used.  This  value  is  good  for  iso¬ 
tropic  materials  and  transversely  isotropic  materials.  Therefore,  the 
program  is  limited  in  that  it  cannot  have  a  general  orthotropic 
material.  However,  the  types  of  materials  that  are  transversely  iso¬ 
tropic  includes  fiber-matrix  materials  and  many  others. 

Implementation  of  Elastic-Plastic  Analysis  in  the  Finite  Element 
Program 


In  incorporating  the  theory  for  plastic  flow  into  the  existing 
elastic  finite  element  program  two  simplifications  were  used. 

The  first  simplification  was  assuming  the  uniaxial  behavior  of 
the  material  was  bilinear.  Figure  14  shows  a  bilinear  material.  By 
doing  this  and  c  are  reduced  to  constants  independent  of  the  stress 
state.  This  obviously  leads  to  much  simplification  in  the  programming. 
The  bilinear  model  is  fairly  accurate  for  most  materials  if  the  plastic 
deformation  is  not  too  large. 

The  second  simplification  is  that  the  stresses  were  not  forced  to 
be  on  the  yield  surface  during  an  entire  load  increment.  Thus,  the 
consistency  equation  given  in  equation  (4.10)  is  violated.  At  the 
start  of  a  load  step  an  element  is  either  plastic  or  elastic.  If  it 
is  elastic,  it  is  considered  to  be  elastic  for  the  entire  load  step. 

If  it  is  plastic,  all  terms  in  equation  (4.18)  are  computed  using 
values  of  O—  at  the  start  of  the  load  step.  At  the  end  of  the  load 
step  the  stresses  are  checked  for  each  element  to  see  if  it  has  yield¬ 
ed.  The  element  is  treated  as  either  elastic  or  plastic  in  the  next 
load  step  according  to  this  check.  Thus,  the  stresses  may  overshoot 
the  actual  yield  surface  for  a  given  time  step.  This  effect  can  be 
minimized  if  load  steps  are  kept  small.  The  alternative  to  this  would 
be  to  do  a  solution  based  on  the  initial  stresses,  check  the  yield 
condition  at  the  end  and  iterate  to  make  sure  the  stress  remains  on 
the  yield  surface.  This  would  cause  a  substantial  increase  in  the 
execution  time  of  the  program.  Since  this  analysis  is  by  its  very 
nature  approximate,  the  substantial  increase  in  execution  time  needed 
to  obtain  the  more  exact  answers  by  iteration  seems  impractical. 
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V.  NUMERICAL  RESULTS  OF  ELASTIC -PLASTIC  ANALYSIS 


Axisymmetric  Perfectly  Plastic  Problem 

The  first  test  of  the  elastic-plastic  numerical  analysis  is  that 
the  numerical  solution  should  match  an  analytic  solution  for  an  axi¬ 
symmetric  problem. 

By  letting  E^  =  .001  E  in  equation  (4.35)  a  perfectly  plastic 
material  was  approximated.  Using  the  example  of  an  isotropic  axisym¬ 
metric  disk  with  a  hole  in  the  center  under  internal  pressure,  a 
numerical  solution  was  obtained.  This  solution  was  compared  to  an 
analytic  solution  of  the  same  problem^0).  There  is  one  difference 
between  these  solutions;  namely,  the  numerical  solution  uses  a  von 
Mises  yield  condition  while  the  analytic  solution  uses  a  Tresca  yield 
condition. 

Allowing  for  the  difference  in  yield  conditions  the  numerical 
solution  agrees  well  with  the  analytic  solution.  In  Figure  15  it  can 
be  seen  that  the  elastic-plastic  boundary  for  both  solutions  has  the 
same  shape  and  that  yielding  at  a  given  radial  location  occurs  at  a 
higher  pressure  for  the  numerical  solution.  This  is  to  be  expected 
from  the  difference  in  the  yield  conditions  (i.e.,  for  this  example 
the  Tresca  condition  does  predict  yielding  at  a  lower  pressure  than 
the  von  Mises  condition) . 

The  circumferential  strain  at  the  outer  surface  is  shown  a  func¬ 
tion  of  pressure  in  Figure  16.  Again  the  analytic  and  numerical 
solutions  agree  well.  Since  there  is  more  yielding  at  a  given  pres¬ 
sure  with  the  Tresca  condition,  the  structure  is  less  stiff.  Larger 
strains  are  to  be  expected  with  the  analytic  solution  at  a  given 
pressure  than  with  the  numerical  solution.  This  is  shown  by  Figure 
16. 

Nonaxisymmetric  Perfectly  Plastic  Problems 

In  doing  the  numerical  solution  in  the  axisymmetric  example  two 
segments  were  used  in  the  6  direction.  Each  segment  had  a  ratio  of 
outer  radius  to  inner  radius  of  2.  Eight  elements  in  the  radial 
direction  and  one  in  the  axial  direction  were  used. 

Two  nonaxisymmetric  cases  were  also  studied.  Both  cases  used  two 
segments  as  before.  In  each  case  one  segment  was  the  same  as  the  axi¬ 
symmetric  configuration.  One  configuration  removed  the  outermost 
element  from  the  other  segment  leaving  seven  elements  in  the  radial 
direction.  The  second  nonaxisymmetric  configuration  had  only  the  inner 
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four  elements  in  the  radial  direction.  These  configurations  are 
shown  in  Figure  17. 

For  both  these  configurations,  a  numerical  solution  was  obtained. 
Since  the  material  is  perfectly  plastic,  the  solution  fails  when  the 
elastic- plastic  boundary  reaches  the  outer  boundary  of  the  structure. 
The  radial  displacement  of  the  inner  boundary  at  the  interface  be¬ 
tween  the  two  segments  is  shown  in  Figure  18.  From  this  graph  it  can 
be  seen  that  as  material  is  removed  from  one  of  the  segments  the 
structure  loses  stiffness  and  there  are  greater  displacements  at  a 
given  pressure.  While  there  is  no  analytic  solution  for  comparison 
m  the  nonaxisymmetric  case,  these  displacements  do  exhibit  expected 
trends. 


Hardening  Cases,  Nearly  Elastic  Examples 

The  final  example  involved  the  configurations  denoted  as  case  1 
and  case  3  in  Figure  17.  These  examples  were  solved  numerically  with 
CT  -  • 9E  in  equation  (4.20).  For  this  value  of  the  tangent  modulus 
the  numerical  solutions  would  not  be  expected  to  vary  greatly  from 
the  elastic  solutions. 


For  the  axisymmetric  case.  Figure  15  shows  how  the  elastic-plastic 
boundary  grows  much  slower  than  in  the  perfectly  plastic  case.  In 
Figure  16  it  can  be  seen  that  the  circumferential  strain  as  a  function 
of  pressure  is  very  nearly  linear.  This  is  consistent  with  the  solu¬ 
tion  being  close  to  the  linear  elastic  solution. 


In  Figure  19  the  displacements  of  the  inner  surface  at  the  inter¬ 
face  between  the  two  segments  are  shown.  It  can  be  seen  that  these 
displacements  are  also  nearly  linear  as  expected.  The  elastic  solu¬ 
tion  is  shown  for  comparison. 
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VI.  CONCLUSIONS 


The  purpose  of  this  investigation  was  to  develop  a  method  that 
could  be  used  to  analyze  nonaxisymmetric  configurations  without  per¬ 
forming  a  full  3-dimensional  analysis.  The  analysis  was  broken  into 
two  parts:  elastic  and  elastic-plastic. 

The  elastic  analysis*  gives  accurate  results  for  a  number  of  ex¬ 
amples.  These  results  were  verified  by  comparison  to  results  obtained 
from  another  finite  element  code  of  proven  accuracy.  Also  these  re¬ 
sults  are  vastly  superior  to  results  obtained  by  the  substructure 
method.  The  substructure  method  is  another  method  which  has  certain 
geometric  similarities  to  the  perturbation  method  presented  in  this 
paper.  Also  the  versatility  of  the  method  was  demonstrated  on  the  tube 
example.  This  showed  that  more  complex  configurations  can  be  handled 
by  the  perturbation  method. 

For  the  elastic-plastic  version  the  analysis  is  designed  to  handle 
both  isotropic  and  orthotropic  materials.  Since  the  examples  presented 
are  all  trying  to  establish  the  accuracy  of  the  method,  only  problems 
with  some  means  of  judging  the  answers  are  presented.  Since  there  are 
no  analytic,  experimental,  or  numerical  results  available  for  ortho¬ 
tropic  materials,  no  orthotropic  examples  are  presented.  Obtaining 
such  results  experimentally  was  beyond  the  scope  of  this  work. 

However,  for  the  isotropic  elastic-plastic  examples  the  results 
are  accurate.  In  this  case  an  analytic  solution  was  available  for 
comparison.  These  examples  represent  the  wide  variety  of  material 
behavior  that  can  be  approximated.  Examples  range  from  perfectly 
plastic  to  nearly  elastic. 

In  conclusion,  the  method  presented  here  gives  an  accurate, 
efficient  means  of  analyzing  a  wide  class  of  problems.  The  method 
is  quite  versatile  in  the  geometric  configurations  it  can  analyze  as 
well  as  the  material  behavior. 
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Figure  1.  Typical  nonaxisymmetric  configuration  to  whic^  the 
structural  analysis  method  can  be  applied. 
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Figure  4.  Subdivision  of  quadrilateral  elements  into  triangular 
elements  and  definition  of  nodal  displacements. 
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Figure  5.  Relation  between  cylindrical  and  orthotropic  coordinates. 
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Figure  6.  End  view  of  nonaxisymnetric  disk  showing  segment 
and  element  numbering  systems. 
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Radial  Displacement  (in 


Figure  7.  Radial  displacement  for  nonaxisymmetric  disk  obtained 
from  different  methods  of  analysis. 
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Figure  8.  Circumferential  stress  for  nonaxisymmetric  disk 
obtained  from  different  methods  of  analysis. 
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Figure  9.  Radial  stress  for  nonaxisymmetric  disk  obtained 
from  different  methods  of  analysis. 
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Figure  10.  Nonaxi symmetric  tube  configuration. 
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Figure  11.  End  view  of  nonaxisyimnetric  tube  configuration  showing 
element  and  segment  numbering  systems. 
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figure  12.  Location  of  three  sets  of  connecting  nodes 
shown  in  r-z  plane. 
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Figure  13.  Kinematic  or  Prager's  strain  hardening. 
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Figure  14. 


Uniaxial  tensile  stress-strain  curve  for  a 
bilinear  material.-- 
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Figure  15„  Radial  location  of  elastic-plastic  boundary  as  a  function  of  pressure 
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Figure  17.  Configurations  analyzed  shown  in  r-z  plane. 
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Figure  18 .  Radial  displacement  at  inner  boundary  as  a  function  of  pressure 
for  perfectly  plastic  materials. 
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Figure  19.  Radial  displacement  at  inner  boundary  as  a  function  of  pressure 
for  strain-hardening  materials. 


TABLE  1 


Stresses  from  various  methods  of  analysis 
for  axisymmetric  disk  under  internal  pressure  of  1000  psi 

Radial  Stress,  0 

rr 


Element  No. 

Axisymmetric 

Perturbation 

Substructure 

i 

-598 

-610 

-  39 

2 

-242 

-243 

-149 

3 

-  97 

-  94 

-199 

4 

-  24 

-  18 

-222 

Circumferential 

Stress,  0QQ 
yy 

Element  No. 

Axisymmetric 

Perturbation 

Substructure 

i 

857 

853 

908 

2 

491 

487 

527 

3 

345 

345 

327 

4 

271 

274 

209 

53 


TABLE  2 


Circumferential  stress  from  various  methods  of  analysis  of  a 
disk  with  internal  pressure  of  1000  psi 


Element  No. 

Perturbation 

Plane  Stress 

Substructure 

1 

1164 

1216 

1193 

2 

505 

540 

544 

Segment  1 

3 

208 

215 

202 

4 

46 

20 

1 

17 

1078 

1133 

1159 

Segment  6 

18 

906 

864 

829 

54 


TABLE  3 


Radial  stress  a  for  nonaxisymmetric  disk  under 
internal  pressure  of  1000  psi 


Element  No. 

Perturbation 

Plane  Stress 

Substructure 

i 

-446 

-519 

109 

2 

-164 

-150 

-  79 

3 

-  53 

-  30 

-164 

Segment  1 

4 

6 

0 

-203 

17 

-608 

-541 

-239 

18 

-  92 

-116 

-325 

Segment  6 

TABLE  4 


Circumferential  stress  from  perturbation  and 
substructure  analysis  for  nonaxisymmetric  tube 
under  internal  pressure  of  1000  psi 


Element  No. 

Perturbation 

Substructure 

Limit  Solution 

i 

1010 

1061 

1164 

2 

503 

544 

505 

3 

284 

268 

208 

4 

108 

102 

46 

17 

1155 

1217 

1078 

18 

811 

760 

906 

61 

934 

986 

857 

62 

486 

528 

491 

63 

301 

286 

345 

64 

205 

141 

271 

81 

755 

956 

857 

82 

479 

594 

491 

83 

386 

402 

345 

84 

345 

288 

271 

56 


TABLE  5 


Radial  stress  a  from  perturbation  and 
substructure  analysis  of  nonaxisymmetric 
tube  under  internal  pressure  1000  psi 


Element  No. 

Perturbation 

Substructure 

Limit  Solution 

1 

-523 

+  51 

-446 

2 

-200 

-104 

-164 

3 

-  71 

-178 

-  53 

4 

-  6 

-216 

+  6 

17 

-574. 

+102 

-608 

18 

-112 

+  6 

-  92 

61 

-570 

+  5 

-599 

62 

-228 

-130 

-243 

63 

-  88 

-193 

-  97 

64 

-  16 

-223 

-  24 

81 

-659 

-106 

-599 

82 

-266 

-212 

-243 

83 

-  98 

-262 

-  97 

84 

-  17 

-286 

-  29 

57 


1 

2 

3 

4 

17 

18 

61 

62 

63 

64 

81 

82 

83 

84 


TABLE  6 


Circumferential  stress  cfgg  for  nonaxisymmetric 
tube  under  internal  pressure  of  1000  psi  with 
various  choices  of  connecting  nodes 


Set  a 

Set  b 

Set  i 

1010 

1011 

1011 

503 

503 

503 

284 

284 

284 

168 

167 

168 

1155 

1155 

1156 

811 

810 

809 

934 

933 

934 

486 

486 

486 

301 

301 

300 

205 

206 

205 

755 

756 

756 

479 

479 

479 

386 

386 

385 

345 

345 

344 
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APPENDIX  A 


DEVELOPMENT  OF  MATRICES  USED  IN  ANALYSES 


A. 1  Derivation  of  Matrix  [B] 

The  [B]  matrix  introduced  in  equation  (2.3)  is  obtained  by  com¬ 
bining  equations  (2.1)  and  (2.2).  Writing  equation  (2.1)  in  matrix 
form  for  the  nodal  displacements  of  a  triangular  element  gives: 


where  RZj[  is  a  row  matrix  [1,  r^,  z-jj  and  i  and  j  vary  from  1  to  3 
and  from  1  to  9  respectively.  The  parameters  (r^,  z^)  are  the  nodal 
coordinates  of  the  triangular  element.  Equation  (A.l)  can  be  ex¬ 
pressed  as: 


{6}  =  [C]  {a}  (A. 2) 

where  {6}  is  the  nodal  displacement,  {a}  is  the  generalized  coordinates 
and  [C]  is  defined  from  equation  (A.l).  Solving  equation  (A. 2)  for 
{a}  gives: 


{a} 


[C]"1  {6} 


(A.  3) 


Then  substituting  equation  (2.1)  into  equation  (2.2)  gives: 


< 


ai 

e 

r 

'o 

1 

0 

0 

0 

0 

0 

0 

0 

a2 

e 

z 

0 

0 

0 

0 

0 

1 

0 

0 

0 

a3 

£0 

>- 

r 

1 

z 

r 

0 

0 

0 

0 

0 

0 

< 

a4 

a5  } 

Yrz 

0 

0 

1 

0 

1 

0 

0 

0 

0 

a6 

YZ0 

0 

0 

0 

0 

0 

0 

0 

0 

1 

a7 

(  -< 

CD 

l _ 

0 

0 

0 

0 

0 

0 

r 

0 

-z 

r 

00  CTi 

(A. 4) 
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liquation  (A.  4)  can  be  expressed  as: 


{e}  =  [Q]  {a}  (A. 5) 

where  [Q]  is  defined  from  equation  (A. 4).  Substituting  for  {a}  from 
equation  (A. 3)  gives: 

{£}  =  [Q]  [C]"1  {6}  (A. 6) 

Making  the  substitution: 

[B]  =  [Q]  [C]"1  (A. 7) 

gives  the  result: 

{e}  =  [B]  {6}  (A. 8) 

which  is  identical  to  equation  (2.3). 


A. 2  Derivation  of  Matrix  [B^] 


In  forming  the  nonaxisymmetric  stiffness  matrix  the  [62]  matrix 
was  introduced  in  equation  (2.13).  This  matrix  is  found  in  a  manner 
similar  to  the  [B]  matrix  explained  in  the  previous  section.  Equation 
(2.10)  can  be  expanded  as: 


where  matrix  C  is  defined  in  equation  (A. 2),  0  is  a  multiplying  con¬ 
stant,  i  varies  over  integers  1  to  6,  and  j  varies  over  integers  1  to 
18.  We  define  0  =  0  on  one  face  of  the  segment  and  0  =  0  on  the  other 
face.  Equation  (A. 9)  can  be  written  in  partitioned  form  as: 


(A. 10) 
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Ot 

where  [C]  is  identical  to  the  [C]  in  equation  (A. 2)  and  (g)  are  the 
IS  generalized  coordinates.  Solving  for  {^}  gives: 


(A. 11) 


where : 


[cx] 


-1 


Substituting  equation  (2.10)  into  (2.11) 


gives : 


(A. 12) 


0,  1,  8x0,  0,  7x0 
5x0,  1,  8x0,  0,  3x0 

l,Zirn6az6_nl1Z 

r>  j. >  0x0,  r  >  3x0,  1,  ^ 

2x0,  1,  0,  1,  6x0,  0,  0,  0,  4x0 
8x0,  1,  3x0,  p  1,  p  2x0,  0 
6x0,  0,  ^ ,  1,  3x0,  0,  ^ 


(A. 13) 


where  j  varies  over  integers  1  to  18  and  the  rectangular  matrix  has 
been  written  in  a  compressed  form  where  nxO  indicates  n  consecutive 
elements  which  are  identically  equal  to  zero.  Equation  (A. 13)  can  be 
rewritten  as: 

{£>  -  [Qj]  {p  (A- 14) 


Combining  equations  (A. 11)  and  (A. 13)  gives: 


(e)  = 

[Q:]  [c^-1  {6} 

(A. 15) 

or 

(e)  = 

[B:]  (6} 

(A, 16) 
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where 


[V  =  tv  [cx] 


(A. 17) 


A.  3  Description  of  [D]  Matrix 

The  [D]  matrix  is  defined  as: 

{a}  =  [D]  {e}  (A. 18) 

The  inverse  of  equation  (A. 18)  is: 

{£}  =  [D]'1  {a}  (A. 19) 

The  [D]  matrix  is  defined  in  equation  (2.4  ).  It  can  be  expressed 
as : 


[D] 


-1  = 


1_ 

E 


sn 


V 


tn 


ns 


1_ 

E 

c 

V 


ts 

"t 

0 


nt 


st 


0  0  0 


0  0  0 


0  0  0 


0  0 


ns 


G  , 
st 


0  0 


tn 


u 


(A. 20) 


The  [D]  matrix  is  obtained  by  inverting  equation  (A. 20) 


A. 4  Derivation  of  the  Matrices  _LBJ_  and  [G] 

If  we  number  the  connecting  nodes  of  a  segment  from  1  to  8  where 
the  nodes  1-4  are  on  one  face  and  the  nodes  5-8  on  the  other  face, 
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the  [R]  matrix  in  the  equation 


<V  =  [R]  fp} 

can  be  partitioned  as: 


[R]  = 


R,  . 


0 


where  [R^]  has  the  form: 


[Rx]  = 


RR 

rr2 

RR, 

RR, 


where  matrices  RR.  are  defined  as  follows: 

1  _ 


[RR,] 


Rzz. 

1 


Rzz.  Rzz. 

l  l 


where  Rzz^  is  a  row  matrix  [1,  r^,  z^,  r^z^]. 


(A. 21) 


(A. 22) 


(A. 23) 


(A. 24) 


The  matrix  [R2]  has  the  same  form  with  rp  replaced  by  r^+^  and  z^  re¬ 
placed  by  zi+4  (i=l,4).  Also  for  each  element  we  have: 


{up}  =  [G]  {p} 


(A. 25) 


where  [G]  can  be  partitioned  as: 


[G] 


(A. 26) 


In  this  case,  since  the  nodes  on  opposite  faces  of  the  segment 
correspond,  their  sub-matrices  are  identical.  The  matrix  G^  has  the 
same  form  as  R4  in  equation  (A. 23)  except  the  coordinates  apply  to 
the  coordinates  of  the  four  corners  of  the  quadrilateral  element  rather 
than  the  four  connecting  points. 
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A. 5  Derivation  of  the  [A]  Matrix 


The  concentrated  loads  such  as  pressure  loadings  are  input  into 
the  program  as  occurring  between  two  nodes,  i  and  j.  The  perturbation 
displacements  at  these  nodes  are  related  to  the  perturbation  displace¬ 
ments  at  the  connecting  nodes  by: 


{uPi}  =  [Ai  {p}  (A-27) 

Substituting  for  {p}  from  equation  (2.17)  gives: 

{upi}  =  [A]  [R]_1  {uCp}  (A- 28) 

where  {u  represents  the  perturbation  displacements  of  node  "i".  A 
similar  expression  is  then  written  for  node  "j". 


If  {g}  is  a  matrix  of  applied  loads  at  node  "i"  then: 


{fi} 


[A. 29) 


where  gr,  gz  and  gg  are  the  components  of  the  load  in  the  three  direc¬ 
tions.  These  components  are  repeated  to  indicate  that  the  load  is 
applied  equally  to  each  face  of  the  segment. 


The  [A]  matrix  in  equation  (A. 26)  can  be  written  in  partitioned 
form  as : 


[A] 


where  [A*]  is: 


A*  '  0 

_l 

0  '  A* 

I 


(A. 30) 


[A*]  =  [RR.] 


(A.  31) 


where  [RR. ]  was  defined  in  equation  [A. 24)  and  where  coordinates  r. , 
refer  to  points  where  the  load  was  applied.  1 
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Then  balancing  virtual  work  gives: 


but 


6u^  {fc}  =  Y  6  uT.  {g} 
P  Pi 


<5u\  =  ([A]  [R]"1  6  u  )T 

pi  kL  j  L  J  Cp^ 


or 


6uT.  =  6  uT  [R]"1  [A]T 

pi  cp  J  L  J 

Then  equation  (A.  30)  becomes: 

T 

6uP  {fc}  =  6  iPp  Z  [R]”1  [A]T  {g} 


where  the  summation  is  over  all  nodes  with  an  external  load 


(A. 32) 

(A. 33) 

(A.  34) 

(A. 35) 
applied 
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APPENDIX  B 


INPUT  CARDS  FOR  COMPUTER  PROGRAMS 


1.  SEGMENT  CONTROL  CARD  (one  for  each  version) 


Elastic  Version: 

Format  (2110) 
Columns  1-10 

11-20 

Elastic-Plastic  Version: 
Format  (2110) 
Columns  1-10 

11-20 


NTYPS  (Number  of  different  types  of 
segments;  4  maximum) 

NTOTS  (Number  of  total  segments; 

8  maximum) 


NTOTS  (Number  of  segments,  8  maximum) 
NOLINC  (Number  of  load  increments) 


2.  SEGMENT  DATA  CARDS 


One  card  for  each  type  of  segment. 

Format  (F10.5,  110) 

Columns  1-10  THETA  (Angle  subtended  by  segment) 

11-20  NST  (The  number  of  segments  of  each 

type;  5  maximum;  not  needed  in 
elastic-plastic  version) 


3.  SEGMENT  NUMBERING  CARDS  (not  needed  in  elastic-plastic  version) 
One  card  for  each  type  of  segment. 

Format  (5110) 
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Columns  1-10 

NUMS  (1) 

(Reference  numbers  for 
segments  of  each  type  in 
global  numbering  system; 
there  can  be  one  to  five 

41-50 

NUMS  (5) 

numbers  per  card  depending 
on  how  many  segments  there 
are  of  each  type.) 

CONNECTING  NODES  CARDS 

One 

card  for  each  segment. 

These  cards  must  be  in  order  accord- 

ing 

to  the  global  numbering 

system  for  segments. 

Format  (8110) 

•  "  ♦'*  ■-< 

Columns  1-10 

NPC  (1) 

(Nodal  number  for  connecting 
nodes  according  to  the 
axisymmetric  grid  for  the 
segment) . 

71-80 

NPC  (8) 

Cards  5-18  must  be  repeated  for  each  different  type  of  segment. 
These  are  the  control  cards  for  the  axisymmetric  solution. 


5.  TITLE  CARD 

Format  (20A4) 

Columns  1-80  TITLE  (Title  for  particular  case) 


6 .  CONTROL  CARD 


Format  (615,  F5.0,  515) 

Columns  1-5  NNLA  (Number  of  nonlinear  approxima¬ 

tions;  NNLA  =  1  for  this  version  of 
the  program) 

6-10  NUMTC  (Number  of  temperature  cards; 

if  -2,  a  constant  temperature  is 
specified) 


11-15  NUMMAT  (Number  of  different  materials; 

6  maximum 


16-20  NUMPC  (Number  of  boundary  pressure 

cards;  200  maximum) 
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21-25 

NUMSC  (Number  of  boundary  shear 
cards;  200  maximum) 

26-30 

NUMST  (Number  of  boundary  shear 
cards  in  tangential  direction; 

200  maximum) 

31-35 

TREF  (Reference  temperature) 

36-40 

INERT  (This  parameter  decides  if 
inertia  loads  will  be  present, 

INERT  +  0  means  zero  values  of  axial 
acceleration,  and  angular  accelera¬ 
tion  and  velocity  for  each  load 
increment) 

51-55 

INCF  (If  INCF  =  0,  then  surface  loads 
for  each  time  increment  will  be  the 
same  as  for  first  increment) 

56-60 

IPLOT  (Plot  parameter,  1  if  plot  re¬ 
quired) 

7.  MESH  GENERATION  CONTROL  CARD 
Format  (515) 


Columns  1-5 

MAXI  (Maximum  value  of  I  in  mesh; 

25  maximum) 

6-10 

MAXJ  (Maximum  value  of  J  in  mesh; 

100  maximum) 

11-15 

NSEG  (Number  of  line  segment  cards) 

16-20 

NBC  (Number  of  boundary  condition 
cards) 

21-25 

NMTL  (Number  of  material  block  cards) 

8.  LINE  SEGMENT  CARDS 

The  order  of  line  segment  cards  is  immaterial  except  when  plots 
are  requested;  in  this  case,  the  line  segment  cards  must  define  the 
perimeter  of  the  solid  continuously.  The  order  of  line  segment  cards 
defining  internal  straight  lines  is  always  irrelevant. 
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4  Circular  arc  specified  by  1st  and 
2nd  points  at  the  ends  of  the  arc 
with  the  coordinates  of  the  center 
of  the  arc  given  as  the  3rd  point 
(delete  I  and  J  for  3rd  point) 

5  Straight  line  as  boundary  diagonal 
for  which  I  of  ist  point  is  minimum 
for  its  row  and/or  I  of  2nd  point 
is  minimum  for  its  row  (input  only 
1st  and  2nd  points) 

6  Straight  line  as  boundary  diagonal 
for  which  I  of  1st  point  and/or 
2nd  point  is  maximum  for  its  row 
(input  only  1st  and  2nd  points) 

NOTE:  In  specifying  a  circular  arc,  the  points  are  ordered  such  that 

a  counterclockwise  direction  about  the  center  is  obtained  upon 
moving  along  the  boundary. 

9.  BOUNDARY  CONDITION  CARDS 


Each  card  assigns  a  particular  boundary  condition  to  a  block 
of  elements  bounded  by  II,  12,  Jl,  J2.  For  a  line  II  =  12  or  J1  = 
J2 .  For  a  point  II  =  12  and  Jl  =  J2. 

Format  (415,  110,  3F10.0) 

Columns 


1-5 

Minimum  I 

6-10 

Maximum  I 

11-15 

Minimum  J 

16-20 

Maximum  J 

21-30 

Boundary  condition  code 

31-40 

Radial  boundary  condition 

code,  XR 

41-50 

Axial  boundary  condition , 

XZ 

51-60 

Tangential  boundary  condition  XT 
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Format  (3(213,  2F8.3),  15) 


<• 


Columns  1-3 

4-6 
7-14 
15-22 
23-25 
26-28 
29-36 
37-44 
45-47 
48-50 
51-58 
59-66 
67-71 

If  the  number  in  column  71  is 


I  coordinate  of  1st  point 
J  coordinate  of  1st  point 
R  coordinate  of  1st  point 
Z  coordinate  of  1st  point 
I  coordinate  of  2nd  point 
J  coordinate  of  2nd  point 
R  coordinate  of  2nd  point 
Z  coordinate  of  2nd  point 
I  coordinate  of  3rd  point 
J  coordinate  of  3rd  point 
R  coordinate  of  3rd  point 
Z  coordinate  of  3rd  point 
Line  segment  type  parameter 


0  Point  (input  only  1st  point 

1  Straight  line  (input  only  1st  and 
2nd  points) 

2  Straight  line  as  an  internal  diagonal 
(input  only  1st  and  2nd  points) 

3  Circular  arc  specified  by  1st  and 
3rd  points  at  the  ends  of  the  arc 
and  2nd  point  at  the  mid-point  of 
the  arc 


jc 
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If  the  number  in  Columns  21-30  is 


XR 

is 

the 

specified 

R-load  and 

0 

xz 

is 

the 

specified 

Z-load  and 

XT 

is 

the 

specified 

0-load 

XR 

is 

the 

specified 

R-displacement 

and 

1 

XZ 

is 

the 

specified 

Z-load  and 

XT 

is 

the 

specified 

0-load 

XR 

is 

the 

specified 

R-load  and 

2 

XZ 

is 

the 

specified 

Z-displacement 

and 

XT 

is 

the 

specified 

0-load 

XR 

is 

the 

specified 

R-displacement 

and 

3 

XZ 

is 

the 

specified 

Z-displacement 

and 

XT 

is 

the 

specified 

0-load 

XR 

is 

the 

specified 

R-load  and 

4 

XZ 

is 

the 

specified 

Z-load  and 

XT 

is 

the 

specified 

©-displacement 

XR 

is 

the 

specified 

R-displacement 

and 

5 

XZ 

is 

the 

specified 

Z-load  and 

XT 

is 

the 

specified 

©-displacement 

XR 

is 

the 

specified 

R-load  and 

6 

XZ 

is 

the 

specified 

Z-displacement 

and 

XT 

is 

the 

specified 

0-displacement 

XR 

is 

the 

specified 

R-displacement 

and 

7 

XZ 

is 

the 

specified 

Z-displacement 

and 

XT 

is 

the 

specified 

©-displacement 
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NOTE:  All  loads  are  considered  to  be  total  forces  acting  on  one 
radian  segment. 

10.  MATERIAL  BLOCK  ASSIGNMENT  CARD 

Each  card  assigns  a  material  definition  number  to  a  block  of 
elements  defined  by  the  I,  J  coordinates. 

Format  (515,  2F10.0,  215) 


Columns  1-5 

Material  definition  number 
(1  through  6) 

6-10 

Minimum  I 

11-15 

Maximum  I 

16-20 

Minimum  J 

21-25 

Maximum  J 

26-35 

Material  principal  property  inclina¬ 
tion  angle  BETA  which  defines  N-S 
plane  orientation  relative  to  z 
direction  (see  Figure  4) 

36-45 

Material  principal  property  inclina¬ 
tion  angle  ALPHA  which  defines  the 
orientation  of  N-T  plane  relative 
to  r-z  plane  (see  Figure  4) 

46-50 

IANG  (If  IANG  =  0,  then  ALPHA  is  same 
for  total  material  block.  If 
IANG  =  1,  the  ALPHA  varies  in 
sign  in  the  I  direction  from 
element  to  element  every  NANG 
elements.  This  will  allow  for 
equal  but  opposite  helical 
angles . ) 

51-55 

NANG  (Number  of  elements  in  the  I 

direction  with  the  same  ALPHA). 

11.  PLOT  TITLE  CARD* 


Format  (20A4) 

Columns  1-80 

Title  (Title  printed  under  each  plot) 
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12. 


PLOT  GENERATION  INFORMATION  CARD* 


Format  (2F10.U) 

Columns  1-10  RMAX  (Maximum  r  coordinate  of  mesh) 

11-20  ZMAX  (Maximum  z  coordinate  of  mesh) 

*N0TE :  Use  only  if  IPLOT  =  1  (plot  required) 

13.  TEMPERATURE  FIELD  INFORMATION  CARDS 

If  NUMTC  in  columns  6-10  of  the  CONTROL  CARD  is  greater 
than  1,  the  temperature  field  is  given  on  cards.  One  card  must  be 
supplied  for  each  point  for  which  a  temperature  is  specified. 


Format 

(3F10.0) 

Columns 

1-10 

R  coordinate 

11-20 

z  coordinate 

21-30 

Temperature 

If  NUMTC  in  columns  6-10  of  the  CONTROL  CARD  is  -2,  a  constant  tem¬ 
perature  field  is  specified;  the  value  is  given  on  a  single  card. 
Format  (F10.0) 

Columns  1-10  Temperature 

14.  MATERIAL  PROPERTY  INFORMATION  CARDS 

The  following  group  of  cards  must  be  specified  for  each 
material  (maximum  of  6). 

a.  MATERIAL  IDENTIFICATION  CARD 
Format  (215,  2F10.0) 

Columns  1-5  Material  Identification  number 

6-10  Number  of  temperatures  for  which 
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properties  are  given  (12  maximum] 


11-20  Mass  density  of  material  (if 
required) 

21-30  Thermal  expansion  parameter  (If  1, 
free  thermal  expansions  on  the 
material  property  cards;  otherwise, 
coefficients  of  thermal  expansion 
are  on  the  material  property  cards.) 

b.  MATERIAL  PROPERTY  CARDS 

Two  cards  are  required  for  each  temperature. 

First  Card 


Format 

(7F10.0) 

Columns 

1-10 

Temperature 

11-20 

Modulus  of  elasticity 

21-30 

Modulus  of  elasticity 

31-40 

Modulus  of  elasticity 

41-50 

Poisson's  ratio.  vXTO 
*  NS 

51-60 

Poisson's  ratio.  v> 

NT 

61-70 

Poisson's  ratio,  vCT 

Second 

Card 

Format 

(6F10.0) 

Columns 

1-10 

Shear  Modulus ,  GM 

IN  o 

11-20 

Shear  Modulus,  G„ 

O  I 

21-30 

Shear  Modulus , 

31-40 

V  or  “n 

41-50 

“sT  or  “s 

51-60 

aTT  or  aT 

N 
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\ 1  ELD  S1RESS  CARDS  (not  needed  in  elastic  version) 
Format  (7F10.0) 


Columns 


1-10  Yield  stress  in  tension  in  N  direc¬ 
tion 

11-20  Yield  stress  in  tension  in  S  direc¬ 
tion 

21-30  Yield  stress  in  tension  in  T  direc¬ 
tion 

31-40  Yield  stress  in  shear  in  NS  direction 

41-50  Yield  stress  in  shear  in  NT  direction 

51-60  Yield  stress  in  shear  in  TS  direction 

61-70  Hardening  parameter  -  C 


16.  INERTIA  LOAD  CARD 

Format  (3F10.0) 

Starting  with  this  input  card  and  including  the  boundary 
force  cards,  this  data  is  to  be  inputted  as  a  block  for 
each  load  step,  that  is  NLINC  times.  There  are  the  follow¬ 
ing  exceptions  to  this : 

a)  If  INERT  =  0,  then  this  card  is  to  be  omitted 
completely  (no  inertia  load). 

b)  If  INCI  =  0,  then  this  card  is  not  repeated,  but 
appears  in  first  block  only  (the  inertia  loads 
are  constant  for  each  load  step). 

c)  If  INCF  =  0,  then  the  following  boundary  pressure 
and  shear  cards  are  to  be  given  only  for  the  first 
block  and  not  repeated  again  (the  pressure  and 
shear  loads  are  constant  for  each  load  increment) , 
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Columns 


1-10 


ACELZ  (axial  acceleration) 


11-20  ANGVEL  (angular  velocity) 

21-30  ANGACC  (angular  acceleration) 

17.  BOUNDARY  PRESSURE  CARDS 

One  card  is  required  for  each  boundary  element  which  is 
subjected  to  a  normal  pressure,  that  is  the  number  of  these 
cards  is  NUMPC  for  each  load  increment. 

Format  (215,  F10.0) 


Columns 

1-5 

Nodal  point  M 

6-10 

Nodal  point  N 

11-20 

Normal  pressure 

As  shown  in  the  figure  below,  the  boundary  element  must  be 
on  the  left  when  progressing  from  M  to  N.  Surface  normal 
tension  is  input  as  a  negative  pressure. 
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18.  BOUNDARY  SHEAR  CARDS 


One  card  is  required  for  each  boundary  element  which  is 
subjected  to  surface  shear,  that  is,  the  number  of  these 
cards  is  NUMSC  for  each  load  increment. 

Format  (215,  F10.0) 

Columns  1-5  Nodal  point  M 

6-10  Nodal  point  N 
11-20  Surface  shear 

As  shown  in  the  figure  below,  the  boundary  element  must 
be  on  the  left  when  progressing  from  M  to  N.  The  positive 
sense  of  the  shear  is  from  M  to  N. 
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19.  BOUNDARY  TRANSVERSE  SHEAR  CARDS 


One  card  is  required  for  each  boundary  element  which  is 
subject  to  transverse  shear,  that  is,  the  number  of  these 
cards  is  NUMSC  for  each  load  increment. 

Format  (215,  F10.0) 


Columns 

1-5 

Nodal  point  M 

X 

6-10 

Nodal  point  N 

11-20 

Surface  transverse  shear 
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20.  BOUNDARY  CONDITION  CONTROL  CARD 


Format  (15) 

Columns  1-5  NRDF  (This  parameter  is  equal  to  the  total 

number  of  displacement  components 
specified  at  connecting  points,  for 
example,  if  one  displacement  component 
was  specified  at  each  connecting  point 
then  NRDF  would  be  equal  to  the  number 
of  connecting  points) 


21.  BOUNDARY  CONDITION  CARDS 

There  are  NRDF  of  these  cards. 

Format  (110,  F10.0) 

Columns  1-10  NREQ  (The  location  of  the  equation  to  be 

modified  in  the  assembled  matrix 
relative  to  the  connecting  points. 

For  example,  if  there  are  12  connecting 
points  and  at  the  fifth  connecting 
point  the  second  component  of  displace¬ 
ment  is  specified  then  this  integer 
would  be  equal  to  3  x  (5-1)  +  2  =  14) 
11-20  U  (The  actual  boundary  condition  valve 

to  be  specified  in  position  NREQ  in  the 
matrix  equations) 
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APPENDIX  C 


PROGRAM  LISTING 

PROGRAM  NONAXI ( INPUT , OUTPUT , TAPE5=±NPJT , TAPE6=OUTPUT , TAPE1 , 

1  TAPE2 , TAPEi , TAPE22 , TAPE23 , TAPE24 , 

2TAPE21 , TAPE25 , TAPE26 ) 

INTEGER  CODE 

COMMON /NPDAT A/ K( 1000) , CODEC  1000) , XR( 1000) ,  Z( 1000 ) , XZ( 1000 ) , 

1NPNUM( 25,80) ,T(1000) ,XT(1000) 

COMMON /ARG/RRR (5) ,ZZZ(5) ,RR(4) ,ZZ(4) ,S(15, 15) , P ( 1 5 ) ,TT(6) , 

1H(6 , 15) , CRZ(6 , 6) ,XI( 10 ) , ANGLE(4) ,SIG( 1 8 ) , EPS( 1 8) , N 
C0MM0N/ELDATa/BETA(1000) ,EPR(1000) , PR (200) ,SH(200) ,1X0000,5) , 

IIP (200) , JP(200) , IS (200) , JS(200) ,ALPHA( 1000) , 1T(200) , JT( 200) , 

2ST (200 ) 

COMMON / BASIC/ ACELZ , ANGVEL , ANGACC , TREF , VOL , NUMNP , NUMEL , NUMPC , NUMSC , 
1NUMST 

COMMON /NXMESH/THETAN (4) ,NST(4) ,NUM3(4,5) ,NPC(8,8) 

COMMON/ANSI /NJMELS (4) ,NUMNPS(4) 

COMMON/NXDATA/NTP, NTYPS , NTS , NTOTS , GTS 1G (24 , 24 , 4 ) 

COMMON /NONAXi/S 1 ( 30 , 30 ) , PI ( 30 ) , THETA , BS 1 ( 6 , 30 ) 

COMMON /SOLVE/X( 4428) , Y( 4428 ) ,TEM( 4428) , NUMTC ,MBAND 

COMMON /TD/IMIN (100), IMAX (100), JMIN ( 25 ) , JMaX ( 25 ) , MAXI , MAXJ , NMTL , NBC 

COMMON / CON VRG/ 1D0NE 

COMMON /PLANE/ NPP 

COMMON/ HESULT/B3(6, 15),D(6,6),C(6,6) ,AR,BB(6,9) ,CNS(6,6) 
C0MM0N/MATP/R0(6) ,E(12, 16,6) ,EE(l6) ,A0FTS(6) 

DIMENSION  TITLE (20) 

Q#  *********************************** 

C  READ  AND  WRITE  CONTROL  INFORMATION 

C*  *********************************** 

READ(5 , 3000 )  NTYPS, NTOTS 
DO  150  1=1, NTYPS 

150  READ (5 , 3001 )  THETAN ( I ) , NST ( I ) 

DO  151  1=1, NTYPS 
J2  =NST (I ) 

151  READ(5 , 3000 )  (NUMS(I, J) , J=1 , J2) 

DO  152  1=  1, NTOTS 

152  READ(5,3002)  (NPC(I, J) , J  = 1 ,8) 

3000  FORMAT (8X10) 

3001  FORMAT (F10.5,I10) 

3002  FORMAT (8110) 

DO  60  1=21,26 
60  REWIND  I 

WRITE(6,3010) 

3010  FORMATC'I"," SEGMENT  DATA  FOR  NONAXISYMMETRIC  PROBLEM") 

WRITE (6, 301 1 )  NTYPS, NTOTS 

3011  FORMAT  ("  NUMBER  OF  TYPES  OF  SEGMENTS  =  ",I5,//, 

1  "  NUMBER  OF  TOTAL  SEGMENTS  =",I5) 

DO  153  1=1, NTYPS 
WRITE(6 , 3012)  I , THETAN ( I ) , NST ( I ) 

3012  FORMAT ("  ",///,"  SEGMENT  TYPE  =  ",I5/,"  THETA  =  ",F10.5/, 

1  "  NUMBER  OF  SEGMENTS  OF  THIS  TYPE  =  ",I5) 

J2  =  N3T(I) 
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o  o  o  o  o 


WHITE (6 , 30l5)(NUMS(l,J) ,J=1 ,J2) 

301 3  FOKMATC  SEGMENT  NUMBERS  IN  GLOBAL  SYSTEM  ARE  ",5I5) 

15S  CONTINUE 

DO  154  1= 1 , NT0T3 

154  WRITE(6,3014)I, (NPC(I, J) ,J=1 ,8) 

3014  FORMaTC  " , "CONNECTING  NODES  FOR  SEGMENT" , 15 , "  ARE",8I5) 

DO  950  NTP  =  1 , NTYPS 

THETA=  THETAN(NTP)  /57  -  295780 

50  READ (5,1000  ) TITLE , NNLA , NUMTC , NUMMAT , NUMPC , NUMSC , NUMST , TREF 
1 , INERT , NLINC , INCI , INCF , IPLOT 

WRITE ( 6 , 2000 ) TITLE , NNLA , NUMTC , NUMMAT , NUMPC , NUMSC , NUMST , TREF , INERT , 
1NLINC 
NPP=0 

************************************ 

GENERATE  FINITE  ELEMENT  MESH 

************************************ 

100  CALL  MESH 

NUMELS(NTP)  =  NUMEL 
NUMNPS(NTP)  =  NUMNP 
lF  (IPLOT.EQ.1)  CALL  MPLOT 

************************************ 

READ  AND  WRITE  T5MPERATURE  DATA 

C*  *********************************** 

103  IF(NUMTC.EQ.O)  GO  TO  440 

IF.( NUMTC. GT.O)  READ(5,1001)  (X(I) , Y(I) , TEM(I) ,  1=  1 , NUMTC) 
IF(NUMTC.EQ.-2)  CALL  TEM2(NUMNP) 

IF ( NUMTC . EQ . -2 )  GO  TO  440 
MPRINT=0 

DO  210  1=1 , NUMTC 
IF(MPRINT.NE.O)  GO  TO  200 
WRITE(6 , 2001 ) 

MPR1NT=59 

200  MPRINT=MPRINT-1 

210  WRITE(6 , 2002)  X(I) , Y(I) ,TEM(I) 

MPRINT=0 

DO  230  N=1, NUMNP 

IF ( MPRINT . NE . 0 )  GO  TO  220 

WRITE(6 , 2003 ) 

MPRlNT=59 

220  MPRINT=MPRINT-1 

CALL  TEMP(R(N) ,Z(N) ,T(N)) 

230  WRiTE(6 , 2004 )  N  ,R(N) ,Z(N) ,T(N) 

440  MPRINT=0 

DO  460  N=1, NUMEL 

IF ( MPRINT. NE.O)  GO  TO  450 

WRITE(6,2008) 

MPR1NT=59 

450  MPR±NT=MPRINT-1 
II=IX(N , 1 ) 

JJ=IX(N ,2) 
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KK=IX(N,3) 

LL=IX(N ,4) 

C 

C  TEM  iS  TEMPORARX  STORAGE  FOR  ELEMENT  TEMPERATURES 
C 

TEM(N)= (T(II)+T (JJ)+T(KK)+T(LL) )/4 .00 
460  wRXTE ( 6 , 2009 )  N , (IX(N , I) , 1=1 , 5) ,BETA(N) , ALPHA(N) ,TEM(N) 

DO  470  K= 1 , NUMEL 
470  T(K)=TEM(K) 

C*  *********************************** 

C  READ  AND  WRXTE  MATERIAL  PROPERTIES 

C*  ********##########################£ 

500  CONTINUE 

DO  510  M=1,NUMMAT 

READ(5, 1004)  MTXPE, (NT, RO(MTXPE) , AOFTS(MTXPE) ) 

WHITE ( 6 , 20 1 0 )  MT  X PE , NT , RO ( MT XPE ) 

READ(5 , 1005) ( (E(X, J,MTXPE) , J=1 , 1 4 ) , 1= 1 , NT) 

IF(AOFTS(MTXPE) .NE. 1 . )  WRITE (6 , 201 1 ) ( (E(I , J ,MTXPE) ,  J=1 , 1 3) ,1=1  NT) 
IF ( AOFTS (MTXPE ) .EQ. 1 , )  WRITE (6 , 201 2) ( (E(I , J ,MTXPE) , J= 1 , 1 j) ,1=1 , NT) 
DO  510  I=NT ,  12 
DO  510  J= 1 , 1 6 

510  E(I, J, MTXPE) =E(NT,J, MTXPE) 

DO  900  NL=1,NL±NC 
WR±TE(6,20i0)  NL 
ACELZ=0.00 

ANGVEL=0 . 00 
ANGACC=0 . 00 

XF(INERT  .EQ.  0)  GO  TO  511 

IF(NL  .NE.  1  .AND.  INCI  .EQ.  0)  GO  TO  511 

C**************************************************mmmmmmm 

C  REAU  AND  WRITE  DXNAMIC  FORCES 

CVVWttKVJHHilHHHHHHIggggjIjIggjIggggjIggggggggggggggggg  ****#*******»*»######*## 

READ(5,1050)  acelz,  angvel,  angacc 
WRXTE(6 , 2051 )  ACELZ,  ANGVEL,  ANGACC 

511  CONTINUE 

C*  *******#*##***#*##########)(#£££££££ 

C  READ  AND  WRITE  PRESSURE  AND  SHEAR  BOUNDARX  CONDITIONS 

0******************,**,,,,,,,^^^^ 

IF(NL  .NE.  1  .AND.  INCF  .EQ.  0)  GO  TO  700 
600  IF(NUMPC.EQ.O)  GO  TO  630 
MPRXNT=0 

DO  620  L=1 ,NUMPC 

IF ( MPR1NT . NE . 0 )  GO  TO  610 

WRITE (6 ,2013) 

MPRINT=58 

610  MPRINT=MPRXNT-1 

READ(5 , 1006 )  IP(L),JP(L),PR(L) 

620  WR1TE(6 ,2014)  IP(L) , JP(L) ,PR(L) 

630  IF(NUMSC.EQ.O)  GO  TO  701 
MPRINT=0 
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DO  650  L= 1 , NUMSC 
IF(MPRINT.NE.O)  GO  TO  640 
WRITE (6 ,2015) 

MPH1NT=58 

640  MPR±NT=MPRINT-1 

READ(5 , 1006)  I3(L) ,JS(L) ,Sh(L) 

650  WRITE(6 ,2014)  IS(L) , JS(L) ,SH(L) 
ro 1  IF(NUMST . EQ. 0 )  GO  TO  700 
MPRINT=0 

DO  680  L= 1 , NUMST 

IF ( MPRINT . NE . 0 )  GO  TO  670 

WRITE (6, 2025) 

MPRINT=58 

670  MPRINT =MPHir)T- 1 

READ(5 , 1006)  IT (L) ,JT(L) ,ST(L) 

680  WRITE (6 , 2014 )IT(L) ,JT(L) ,3T(L) 

Q*  a****#######**####****#***#*###**** 

C  DETERMINE  BANDWIDTH,  INITIALIZE  ELASTIC-PLamTIC  RATIO, 

C  AND  CONVERT  BETA  FROM  DEGREES  TO  RADIANS 

C*  *********************************** 

?00  J=0 

DO  710  N  = 1 , NUMEL 
IX (N  ,5)=IABS(IX(N , 5) ) 

DO  710  1=1,4 
DO  710  L= 1 , 4 
KK=IABS(IX(N ,I)-IX(N,L)) 

IF(KK.GE.J)  J=KK 
HO  CONTINUE 
MBAND=3*J+3 
IF(NL.GT.I)  GO  TO  721 
DO  720  N=1, NUMEL 
EPR(N)=1 . 

ALPHA(N)=ALPHA(N)/57.2y5780 
720  BETA(N)=BETA(N)/57. 295780 
CONTINUE 

C*  *********************************** 

C  SOLVE  NONLINEAR  PROBLEM  BY  SUCCESSIVE  APPROXIMATIONS 

C*  a********************************** 

DO  800  NNN= 1 , NNLA 

C 

C  FORM  STIFFNESS  MATRIX 

•C 

CALL  STIFF 

SOLVE  FOR  DISPLACEMENTS 
CALL  SOLV 
COMPUTE  STRESSES 
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799 

800 
810 

900 

950 


910 

1000 

1001 

1004 

1005 

1006 
1030 
2000 


WRITE (6, 20 16)  NITER 
WRITE(6 , 201 7)  NITER 


2001 

2002 

2003 

2004 
2008 


CALL  STRESS 
CALL  STORE 

IFdDONE . NE . 1 )  GO  TO  800 
NiTER=NNN 

IF(IOONE.EQ. 1 )  GO  TO  810 
CONTINUE 
IFdDONE.EQ.  1 ) 

IF(IDONE.NE.I) 

CONTINUE 
CONTINUE 
CALL  ASEMbL 
CALL  ANSWER 
CONTINUE 

FOKMAT(20A4/6I5,F5. 0,515) 

FORMAT(3F10 . 0) 

FORMAT  ( 215, 2F 10.0) 

FORMAT (7F 10.0) 

FORMAT  (2I5.F10.0) 

FORMAT(SFIO.O) 

FORMAT  (2H1  ,20A4/ 

NUMBER  OF  APPROXIMATIONS - 14/ 

TEMPERATURE  CARDS— 14/ 

MATERIALS - 14/ 

PRESSURE  CARDS - 14/ 

SHEAR  CARDS - 14/ 

TORSION  CARDS - 14/ 

TEMPERATURE - El  2. 4/ 

INERTIA  CARDS - 14/ 

LOAD  INCREMENTS - 14/) 

(1H1 , 1 jX, 1HR, 14X, 1HZ, 14X, 1HT) 

(3F15.3) 

(35H1  NR  -  Z 

( 15 , 2F10 . 4 ,F10 . 0) 

(74H1  EL  I  J 
TEMPERATURE) 

(15,414, 18, F1 1 . 1 ,2F13.3) 

(1  HI ."MATERIAL  IDENTIFICATION  NUMBER  =”,12/ 
, "NO .  OF  MATERIAL  TEMPERATURE  CARDS  =",12/ 

, "MASS  DENSITY  =”,E15.7) 


33H0 
33H0 
33HO 
33HO 
33H0 
33H0 
33H0 
33HO 
33H0 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
1LPHA 

2009  FORMAT 

2010  FORMAT 
1 1  ri 
21H 


NUMBER 
NUMBER 
NUMBER 
NUMBER 
NUMBER 
REFERENCE 
NUMBER  OF 
NUMBER  OF 


OF 

OF 

OF 

OF 

OF 


K  L 


•  T) 

MATERIAL  ANGLE  BETA 


2011  FORMAT  (1H  , "TEMPERATURE  =",E15.7/ 

1 1 H  , "MODULUS  OF  ELASTICITY-EN  =”,E15.7/ 

21 H  ."MODULUS  OF  ELASTIC ITY-ES  =",E15.7/ 

31H  , "MODULUS  OF  ELASTICITY-ET  =",E15.7/ 

41H  , "POISSON  RATIO-NUNS  =”,E15.7/ 

51H  ."POISSON  RATIO-NUNT  =”,E15.7/ 

61H  ."POISSON  RATIO-NUST  =",E15.7/ 

71H  ."SHEAR  MODULUS-GNS  =",E15.7/ 

81H  ."SHEAR  MODULUS-GST  =",E15.7/ 

9 1 H  ."SHEAR  MODULUS-GTN  =”,E15.7/ 

1 1 H  , "COEFFICIENT  OF  THERMAL  EXPANSION-AN  =”,£15-7/ 


ANGLE  A 
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21H  ."COEFFICIENT  OF  THERMAL  EXPANSlON-AS  =",E15.7/ 

J1H  , "COEFFICIENT  OF  THERMAL  EXPANSION-AT  =",E15.7/) 

2012  FORMAT  (1H  , "TEMPERATURE  =",E15.7/ 

1 1 H  ."MODULUS  OF  ELASTICITY-EN  =",E15.7/ 

21H  ."MODULUS  OF  ELASTICITY-ES  =",E15.7/ 
ilH  ."MODULUS  OF  ELASTICITY-ET  =",E15.7/ 

41H  ."POISSON  RATIO-NUNS  =",E15.7/ 

51H  , "POiSSON  RATIO-NUNT  =",E15.7/ 

61H  ."POISSON  RATIO-NUST  =",E15.7/ 

7 1 H  ."SHEAR  MODULUS-GNS  =",E15.7/ 

81H  ."SHEAR  MODULUS-GST  =",E15.7/ 

9 1 H  ."SHEAR  MODULUS-GTN  =",E15.7/ 

1 1 H  ."FREE  THERMAL  STRAIN-FN  =",E15.7/ 

21 H  ."FREE  THERMAL  STRAIN-FS  =",E15.7/ 

3 1 H  ."FREE  THERMAL  STRAIN-FT  =",E15.7/) 

20 1^  FORMAT  ( J0H1  PRESSURE  BOUNDARY  CONDITiONS/20H  I  J  PRESSURE) 

2014  FORMAT  (2I5.F10.1) 

2015  FORMAT  (27H1  SHEAR  BOUNDARY  CONDITIONS/ 1 7H  I  J  oHEAR) 

2016  FORMAT  (26H  THE  SYSTEM  CONVERGED  IN  I2.11H  ITERATIONS) 

201 J  FORMAT  (j53H  THE  SYSTEM  DID  NOT  CONVERGE  IN  I2.11H  ITERATIONS) 

2024  FORMAT  (43HO  THE  AXISYMM'ETRIC  OPTION  NAS  BEEN  SELECTED) 

2025  FORMAT(30H1  TORSION  BOUNDARY  CONDITIONS/ 1 7H  1  J  SHEAR) 

2030  F0fiMAT(1H1 ."LOAD  STEP=",I4) 

2031  FORMAT(1riO  ."AXIAL  ACCELERATION  =",E12.4/ 

1 1  HO  ."ANGULAR  VELOCITY  =",E12.4/ 

21  HO  ."ANGULAR  ACCELERATIONS" .E12.4) 

920  STOP 
END 

SUBROUTINE  ANGLE  (R , Z , RC ,ZC , ANG) 

C  FIND  ANGLE  OF  INCLINATION  BETWEEN  0  AND  2*Px 

Q*****#************<##>#>>Jt>)tj(j(JtJtJtJtJt> 

PI=3. 1415927 
D1=(Z-ZC) 

D2=(R-RC) 

IF(ABS(fl-RC) .GT.1.E-8)  GO  TO  100 
ANG=PI/2 . 

IF(D1.GT.1 .E-8)  RETURN 

ANG=-ANG 

RETURN 

C*************#**##*##,###*#*****,,** 

C  ALLOW  CIRCLE  TO  CROSS  AXIS 

(j******************#####Jtl(jtJtl(Jtl(JtJt)(JtJtJt 

100  ANG  =AT AN2 ( D 1 , D2 ) 

RETURN 

END 

SUBROUTINE  ANSWER 
INTEGER  CODE 

COMMON/ ELDATA/BETA( 1000) , EPR( 1 000 ) , PR(200) ,SH(200) ,IX( 1000, 5) , 

IIP (200) ,JP(200) ,IS(200) , JS(200) .ALPHA(IOOO) , IT (200)  ,JT(200)  , 

2ST (200 ) 
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C0MM0N/ARG/RRR(5) ,ZZZ(5) ,RH(4) ,ZZ(4),S(15,15) ,P(15),TT(6), 

1H(6 , 15) ,CHZ(6,6) ,XI( 10) ,ANGLE(4) ,SiG(  1 8 ) , EPS(  1 8) ,N8 
COMMON/NXSOLV/SKG( 1 32 , 24 ) ,FTG( 132), ITOT 
C0MM0N/ANS2/  UT1(24),  G(24,24),  GR 1 (24 , 24 ) , DUMM(24 , 24) 

COMMON / AIMS  1  /NUMELS ( 4 )  ,NUMNPS(4) 

COMMON /NON AXI/81  (30,30)  ,P1 (30) , THETA ,B31 (6,30) 

COMMON /N  XD AT A/NT  P , NTYPS , NTS , NTOTS , GTS 1G ( 24 ,24,4) 
COMMON/NXMESH/IHETAN(4 ) ,NST(4) ,NUMS(4,5) ,NPC(8,8) 

COMMON/AHGI/SiGI (18) ,EPS1 (18) 

COMMON /S0LVE/B( 162) , A( 1 62 , 81 ) ,NUMBLK,MBAND 
DIMENSION  UT(24)  ,UC1 (24) ,UC(24) , HI (24,24) 

REWIND  25 
REWIND  26 
KOLD=1 

DO  100  K=1 , NTYPS 
KNEW=K 
J1=  NST(K) 

NUMNP  =  NUMNPS(K) 

NUMNP3  =  3*NUMNP 
NOMEL  =  NUMELS(K) 

K20  =  K+20 

READ (26 )  (B(I) ,1=1 ,NUMNP3) 

READ(26)  ((IX(I,J) ,J=1 ,4) ,1=1 , NUMEL) 

DO  100  L= 1 , J 1 
NS=NUMS(K,L) 

REWIND  K20 
WRITE(6, 1200)  K,NS 
READ(25)((H1(I,J) , J=1 ,24) ,1=1 ,24) 

DO  110  KK=1 ,4 
NP 1  =  NPC(NS ,KK) 

NP2  =  NPC(NS ,KK+4) 

DO  110  1=1,3 

UC(3*(KK-1 )+I)  =  B(j*NP1-3+I) 

UC(3*(KK-1 )+i+12)  =  B(3*NP2-3+l) 

110  CONTINUE 

DO  115  KK=1 , 24 

115  UT(KK)  =  FTG(KK+(NS-1 ) *1 2 ) 

WRITE(6 , 900 ) 

900  FORMAT ( "  EL  SIGMAR  SIGMAZ  SIGMAC  SIGMARZ  SIGMAZC" 

1  ,"  SIGMACR  SIGEFF",/"  EPSR  ERSZ  EPSC" , 

2  "  EPSRZ  ERPSZC  EPSCR") 

IF(KOLD.EQ.KNEW)  REWIND  21 

IF ( KOLD . NE . KNEW )  KOLD=KNEW 
DO  120  N=1, NUMEL 

READ(K20)((CRZ(I,J) , J=1 ,6) ,1=1 ,6) 

READ(K20) ( (BS1 (I , J) ,J=1 ,30) ,1=1 ,6) 

READ(K20)((  G(I,J) ,J=1 ,24) ,1=1 ,24) 

DO  125  1=1,24 
DO  125  J= 1 , 24 

GR1 (I, J)  =  0.00 
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DO  125  M=1 ,24 

125  OR  1 ( 1 , J )  =  GR1 (I, J)  +  G(I,M)*R1(M,J) 

DO  126  1=1,24 
DCKI)  =0.00 
and)  =o.oo 

DO  126  J=1 ,24 

ucid)  =  ucid)  +  gri  d,  j)*uc(j) 

126  arid)  =  uTid)  +  grki,j)*ut(j) 

DO  130  1=1,4 
l±=3*x 
JJ=3*iX(N,l) 

PKII-2)  =  B(  JJ-2) 

Pldi-1)  =  B(  JJ-1 ) 

P1(II  )  =  B(JJ  ) 

Pldl+10)  =  B(  JJ-2 ) 

PKll+11)  =  B(  JJ-1 ) 

Pldi+12)  =  B(JJ) 

1;j0  CONTINUE 

DO  135  1=1 ,24 

135  P1(I)  =  Pld)  -UC1d)+JT1d) 

DO  136  1=1,3 

P 1(1+24)=  (P 1 (i)+P 1 (1+3 )+P 1 (1+6 )+P 1 (1+9 ) ) /4 . 00 
1 3o  PI (±+2? )  =  (P 1 (1+12)+P1 (1+15 )+P1 (1+13)+P1 (1+21 ) ) /4 .00 
DO  140  1=1,6 
EPS1 (1 )  =  0.00 
DO  140  J=1 ,30 

140  EPSId)  =  EPS1  (D+BS1  (I , J)*P1  (J) 

DO  150  1=1,6 
31G1(1)  =  0.00 
DO  150  J=1 ,6 

150  S1G1(1)  =  SlGI(l)  +  CRZ(l,  J)*EPS1 (J) 

31GEFF=(31G1(1 )-3IG1(2))**2+(SiG1 (2)-SIG1 (3) )**2+(SlG1 (3)- 
1  31G1 ( 1 ) ) **2+6 . * (31G1 (4 )**2+SlG1 (5 )**2+31G1 (6 )**2 ) 

31GEFF=3URT( .5*S1GEFF) 

DO  141  J=1 , 6 

141  EP31  ( J.)  =  EPS1  (J)*100.0 

WR1TE(6 , 1000 )N , (31G1 (1 ) ,1=1 ,6) , S1GEFF 
WRXTE(6 , 1 100) (EPS1 (I)  ,1=1,6) 

120  CONTINdE 
100  CONTINUE 

1000  FORMaT("  ",I5,6F9.0,3X,F9.0) 

1100  FORMATC'  "  ,5X,6F9.5) 

1200  FORMAT ( "  1 " , " SEGMENT  TYPE" ,15,//,"  "/'SEGMENT  NUMBER  =  ",I5) 
RETURN 
END 

SUBROUTINE  A3EMBL 

COMMON/GLB3EG/FI(24,8),FE(24,8),UC(24,8),3K(24,24,8) 
COMMON/NXDATA/NTP , NTYPS , NTS , NTOTS , GTS 1 G ( 24 , 24 , 4 ) 
COi'lMON/NXSOLV/SKG(  1 32,24)  ,FTG(  1  32 ) ,  ITOT 
COMMON/AN32/FC(24) ,G(24,24) , GRI (24,24) ,DUMM(24, 24) 
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ITOT=  24  +  12*(NT0TS-1) 

DO  10  1=1, ITOT 
FTG(i)  =  0.00 
DO  10  J  =  1  ,24 
10  SKG ( 1 , J )  =  0.00 
DO  100  M= 1 , NT0T3 

C  WRITE(6,1000)(FE(I,M)  ,1=1  ,24) 

WRITE(6,1000)(FI(I,M)  ,1=1  ,24) 

WRITE(6,1000)((SK(1,J,M)  ,  J=1  ,24)  ,1=1  ,24) 

WRlTE(6 , 1000 ) ( UC(1 ,M) ,1=1,24) 

1000  FORMAT ("  ",12fi10.3) 

************************************************************** 

COMBINE  FI,  FE,  AND  SK*UC  INTO  A  TOTAL  FORCE  VECTOR  FC 

************************************************************** 

DO  55  1=1,24 
FC(X)  =  0.00 
DO  55  J=1 ,24 

55  FC(i)  =  FC(I)  +  SK(I,J,M)»  UC(J,M) 

DO  60  1=1,24 

60  FC(I)  =  FC(l)  +FE(I ,M)  -FI(I,M) 

if*#######***#####**#****################### «** ************ 

NOW  FILL  GLOBAL  FORCE  AND  STIFFNESS  MATRICES 

C  **** ******************************************* *********** 

DO  70  1=1,24 
II  =  I+(M-1)*12 
FTG(II)  =  FTG(II)  +  FC(1) 

DO  70  J=± , 24 

SKG(I 1 , J+1 -I )  =  SKG  (i1,J+1-I)  +  SK(I , J,M) 

70  CONTINUE 
100  CONTINUE 

C  READ  THE  TOTAL  NUMBER  OF  RESTRAINED  DEGREES  OF  FREEDOM 
READ(5 , 1200)  NRDF 
WHIT£(6 , 1255)  NRDF 

C  J.MPOSE  BOUNDARY  CONDITIONS  ON  RESTRAINED  D-O-F 
DO  150  NBC=1 ,NRDF 

C  READ  THE  EQUATION  NUMBER  AND  THE  IMPOSED  BOUNDARY  CONDITION 
READ (5 , 1250)  NREQ,U 
WRIT£(6, 1260)NREQ,U 
CALL  XMODFY(U,NR£Q) 

150  CONTINUE 
1200  FORMAT (15) 

1250  FORMAT ( 15 , F 1 0 . 0 ) 

1255  FORMAT( 1  HI , "NUMBER  OF  RESTRAINED  DEGREES  OF  FREEDOM  =",110/ 
1  "  EQUATION  NUMBER  VALUE  ") 

1260  FORMAT  (»  " ,5X, I5,5X,F10 .2) 

CALL  XSOLVE 
WRITE(6, 1050) 

WRITE(6 , 1 100 ) ( FTG(l) ,1=1, ITOT ) 

1050  FORMATC'I", "TOTAL  DISPLACEMENTS  AT  CONNECTING  NODES"/ 

1  1 8X , 2HUR , 1 8X , 2HUZ , 1 8X , 2HUT ) 
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1100  FORMAT ( "  ",3E20.7) 

RETURN 

END 

SUBROUTINE  CIRCLE ( ANG 1 , DELPHx , RSTRT , ZSTRT , RC , ZC , I , J ) 

INTEGER  CODE 

COMMON/TD/IMlN (100) ,IMAX( 100) , JMIN(25) ,JMAX(25) , MAXI , MAXJ , NMTL , NBC 
COMMON /NPD AT A/ R( 1 000 ) ,CODE( 1000), XR (1000), Z(1000),XZ  (1000), 

1 NPN  UM( 25 , 80 ) , T( 1 000 ) , XT (1000) 

DIMENSION  AR(25,80) ,AZ(25,80) 

EQUIVALENCE  (R( 1 ) , AR ) , ( Z( 1 ) , AZ) 

C*  a*******##*#**#*#*##### ************ 

C  FIND  INTERSECTION  OF  LINE  AND  CiRCLE  =  NEW  R  AND  Z 

ANG 1 =ANG 1 +DELPHI 

RR=3QRT ( ( RSTRT-RC) **2+( ZSTRT-ZC ) **2 ) 

AR(I, J)=RC+RR*C0S(ANG1 ) 

AZ ( I , J ) =ZC+RR*SIN ( ANG 1 ) 

RETURN 

END 

SUBROUTINE  INTER 

COMMON/ ARG/ RRR(5 ) ,  ZZZ(  5  )  , RR(4 )  , ZZ(4 )  ,S(15, 15)  ,P(15) ,  TT(6 ) , 

1H(6 , 15) , CRZ(6, 6) , XI (10) , ANGLE ( 4 ) , SIG ( 1 8 ) , EPS ( 1 8 ) , N 
COMMON /PLANE/NPP 
DIMENSION  XM(7)»R(7),Z(7)»XX(9) 

DATA  XX/3*. 1259391 805448, 3*. 132394 1527884, .225, 

1  .696140478028, .410426192314/ 

R(7)=(RR(1 )+RR(2)+RR(3) )/3 .0 
Z(7)=(ZZ(1)+ZZ(2)+ZZ(3))/3.0 
DO  100  1=1,3 
J =1+3 

R(I)=XX(8)*RR(I)+(1 . 00-XX(8 ) ) *R(7 ) 
fl(J)=XX(9)*RR(I)+(1 .00-XX(9))*R(7) 

Z(I)=XX(8)*ZZ(I)+(1 . 00-XX( 3 ) ) *Z( 7 ) 

100  Z(J)=XX(9)*ZZ(I)+(1.00-XX(9))*Z(7) 

DO  200  1=1,7 
200  XM( I ) =XX( I ) *R ( I ) 

DO  300  1=1 ,10 
300  XI(I)=0.00 

AREA=.50*(RR(1 )*(ZZ(2)-ZZ(3))+RR(2)*(ZZ(3)-ZZ(1 ) )+Rfl(3)*(ZZ( 1 ) 

1  -ZZ( 2 ) ) ) 

IF(NPP.NE.O)  GO  TO  600 

DO  400  1=1,7 

XI ( 1 ) =XI ( 1 )+XM( I ) 

Xl(2)=XI(2)+XM(I)/R(I) 

XI(3)=XI(3)+XM(I)/(R(I)**2) 

XI(4)=XI(4)+XM(I)*Z(I)/R(I) 

XI(5)=XI(5)+XM(I)*Z(I)/(R(I)**2) 

XI(6)=XI(6)+XM(I)*(Z(I)**2)/(R(I)**2) 

XI (7 )=XI ( l )+XM(I )*R(I) 

XI ( 8 ) = XI ( 8 ) +XM ( I ) *  Z ( I ) 
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XI(9)=X1(9 )+XM(I)#(R(l)**2) 

400  XI(10)=XI(10)+XM(I)*R(I)*Z(I) 

DO  500  1=1 , 10 
500  XI(I)=XI(i) * AREA 
RETURN 

600  XI ( 1 ) =AREA 

Xl(7)=R(7)*AREA 

XI(3)=Z(7)*AREA 

RETURN 

END 

SUBROUTINE  MESH 
INTEGER  CODE 

DIMENSION  AR( 25 , 80 ) , AZ(25 , 80) , NCODE(25 , 80 ) 

COMMON/TD/ IMIN (100), IMAX( 100), JMIN (25 ) , JMAX(25 ) ,MAXI ,MAXJ , NMTL , NBC 
COMMON /NPDAT A/ R( 1000) ,CODE( 1000) ,XR( 1000) ,Z( 1000) ,XZ( 1000) , 
1NPNUM(25,80) ,T( 1000) ,XT( 1000) 

COMMON/ELDATA/BETA( 1000) ,EPR( 1000) , PR (200) ,3H(200) , IX ( 1000,5) , 

IIP (200) , JP(200) , IS (200) , JS(200) , ALPHA ( 1000) , IT (200) , JT(200) , 
2ST(200 ) 

EQUIVALENCE  ( R( 1 ) , AR) , ( Z( 1 ) , AZ) , ( IX( 1 , 1 ) ,NCOD£) 

************************************ 

MESH  CONTROL  INFORMATION 

q*  *********************************** 

READ  (5,1000)  MAXI ,MAXJ,NSEG, NBC, NMTL 
WRITE ( 6 , 2000 )  MAXI , MAXJ , NSEG , NBC , NMTL 

C*  *********************************** 

C  INITIALIZE 

q*  *********************************** 

I3EG=-1 
Pi=5. 1415927 
DO  110  J=1 , 100 
DO  100  1=1 ,25 
NCODE(I, J)=0 
AR(I, J)=0. 

AZ( 1 , J ) =0 . 

JMAX(I)=0 
100  JMIN(I)=MAXi 
IMIN ( J) =MAXJ 
110  IMAX(J)=0 

C*  *********************************** 

C  LINE  SEGMENT  CARDS 

C************************************ 

150  ISEG=ISEG+1 

159  IF ( ISEG . EQ . NSEG )  GO  TO  400 

READ(5 , 1001 )  I1,J1 ,R1,Z1 ,12, J2,R2, Z2, 13, J3,R3,Z3, IPTION 

WRITE(6 ,2001 )I1 ,J1 ,R1 ,Z1 ,12, J2,R2,Z2,Ii,J3,R3,Z3,IPTION 

IPTION=IPTION+1 

AR ( I 1 , J 1 ) =R 1 

AZ(I1 , J1 )=Z1 

NCODE(I1 ,J1 )=1 
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C*  * 

c 

c*  * 

200 


C 


CALL  MNIMX(1 1 , J 1 ) 

GO  TO  (150,200,200,300,300,200,200),  IPTION 

***************,,,,,,,,  #  #  „  #  #  ,  #  #  ,  „  # 

GENERATE  straight  lines  on  boundary 

•  •••a*******#**,,  ***************** 

Di=  ABS ( FLOAT ( 12-1 1 ) ) 

DJ  =  A6S(FL0AT (J2-J1 ) ) 

AR(I2,J2)=R2 
AZ(I2, J2)=Z2 
NCODE(i2, J2)=1 
CALL  MNIMX(T2,J2) 

I3TRT=I1 
ISTP=I2 
JSTRT=J 1 
JSTP=J2 

DIFF=MAX1 (DI,DJ) 

ITEfl=DIFF-1  . 

IINC=0 

JINC=0 

IF(X2.Nfi. 11 )  IINC=(12-1 1 ) /TABS (12-11 ) 

IF(J2.NE. J1 )  JINC= ( J2-J 1 )/IABS(J2-J1 ) 

KAPPA= 1 

iF(I2.NE.I1 .AND. J2.NE. J1 .AND.IPTlON.NE.3)  KAPPA-2 
IF(KAPPA.EQ.2)  DlFF=2.*DlFF 
RINC=(R2-fl1 )/DIFF 
ZiNC=(Z2-Z1)/DlFF 

WRITE ( 6 , 2002 )  D1 , DJ , DIFF , RINC , ZINC , ITER , I INC , JlNC , KAPPA 


CHECK  FOR  INPUT  ERROR 


IF(KAPPA.N£.2.0R.DI.EQ.DJ)  GO  TO  210 
WRiT£(6 , 2003) 

GO  TO  150 


INTERPOLATE 


210  1=11 

J=J1 

WRiTE(6 ,2004) 

DO  230  M= 1 , ITER 

IF (ITER . EQ. 0 . AND . IPTION . EQ. 2 )  GO  TO  230 

IF ( ITER. EQ. 0 .AND . IPT10N . EQ. 6 )  GO  TO  230 

IF ( ITER. EQ.O. AND. IPTION. EQ. 7)  GO  TO  230 

IF (KAPPA. EQ. 2)  GO  TO  220 

IOLD=I 

I=i+IINC 

JOLD=J 

J=J+JINC 

AR ( i , J ) = AR ( IOLD , JOLD ) +RINC 
AZ(I, J)=AZ(IOLD, JOLD)+ZINC 
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WRITE (6, 2005)  1 , J , AH  (I , J) , AZ(i , J) 

CALL  MNlMX(l, J) 

NCODE(I , J)=1 
GO  TO  230 

220  CONTINUE 

lF ( I 1 . GT . 12 . AND . IPTION . EQ . 7 )  GO  TO  221 
IF ( 1 1 . LT . 12 . AND . IPT  ±0N . EQ . 6 )  GO  TO  221 
IOLD=I 
I=I+IINC 

AR(I,J)=AR(IOLD,J)+RINC 
AZ(I,J)=AZ(IOLD,J)+ZINC 
WRITE(6,2005)  I,  J,  Aft  (I,  J)  ,AZ(1 ,  J) 
NCODE(I, J)=1 
CALL  MNIMX(I,J) 

JOLD=J 

J=J+JINC 

AR(l,J)=AR(I, JOLD)+RINC 
AZ(I,J)=AZ(I,JOLD)+ZINC 
NCODE ( I , J ) = 1 

WRITE(6,2005)  i,J,AR(I,J) ,AZ(1, J) 

CALL  MNIMX(I, J) 

GO  TO  230 

221  JOLD=J 
J=J+JINC 

AR (I , J ) = AR (I , JULD ) +RINC 
AZ ( I , J ) =AZ ( I , JOLD ) +ZINC 
NCODE(I,J)=1 

WRITE(6 , 2005 )  1,J,AR(1,J) ,AZ(I, J) 

CALL  MNIMX(I , J) 

IOLD=I 

I=I+I1NC 

AR (I , J ) = AR (TOLD , J ) +RINC 
AZ (I , J ) =AZ ( IOLD , J ) +ZINC 
NCODE(l,J)=1 

WRITE(6 , 2005 )  I, J, Art(If J) ,AZ(I,J) 

CALL  MNIMX(1, J) 

230  CONTINUE 

IF ( KAPPA. EQ. 1 )  GO  TO  150 

IF(I1 .GT. 12. AND. IPTION. EQ. 7)  GO  TO  231 

IFdl.LT. 12. AND. IPTION. EQ. 6)  GO  TO  231 

IOLD=I 

I=I+IINC 

AR (I , J ) =AR (IOLD , J ) +RINC 
AZ( I , J ) =AZ ( IOLD , J ) +ZINC 
GO  TO  232 

231  CONTINUE 
JOLD=J 
J=J+JINC 

AR ( I , J) =AR ( I , JOLD ) +RINC 
AZ(I,J)=AZ(I,JOLD)+ZINC 


93 


o  o 


2^2  CONTINUE 

NCODE(I, J)=1 

WH1TE(6 , 2005 )  I , J , Afi ( I , J ) ,AZ(i,J) 
CALL  MNIMX(I , J ) 

GO  TO  150 

C  GENERATE  CIRCULAR  ARCS  ON  BOUNDARY 

500  AH (12 , JS )=R2 
AZ(I2, J2)=Z2 
NC0DE(I2,J2)  =  1 
CALL  MNIMX(i2,J2) 

IF(IPT10N .EQ. 5)  GO  TO  320 
C 

FIND  CENTER 'OF  CIRCLE 


AR(I3, Ji)=R3 
AZ(I3, J5)=Z5 
NC0DE(I3,J3)=1 
CALL  MNIMX(I3,J3) 

SLAC= ( Z2-Z 1 ) / ( R2-R 1 ) 

SLBFs-1 ./SLAC 
3LCE=(Z3-Z2)/(R3-R2) 

SLDF=-1 . /SLCE 
C 

C  CHECK  FOR  INPUT  ERROR 
C 

IF (ABS (SLAC -SLCE ) .GT . . 001 )  GO  TO  310 
WRITE(6 , 2006 )  R 1 , Z 1 , R2 , Z2 , R3 , Z3 , SLAC , SLCE 
GO  TO  150 

310  R4=R1+(R2-R1 )/2. 

Z4=Z1+(Z2-Z1 )/2 . 

R5=R2+(R3-R2)/2. 

Z5=Z2+(Z3-Z2)/2. 

BBF=Z4-SLBF*R4 
BDF=Z5-SLDF*R5 
RC  = ( BBF -BDF ) / ( SLDF-SLBF ) 

ZC=SLBF*RC+BBF 
WRITE(6 , 2007 )  RC, ZC 
KAPPA= 1 
GO  TO  330 
320  KAPPA=2 
RC=R3 
ZC  =  Z3 

330  ISTRT=I 1 
ISTP=I2 
JSTRT=J1 
JSTP=J2 
RSTRT=R 1 
RSTP=R2 
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ZSTRTsZI 
ZSTP=Z2 

340  CALL  ANGLE(RSTRT,ZSTRT,RC,ZC,ANG1) 

CALL  ANGLE ( RSTP , ZSTP , RC , ZC , ANG2 ) 

IF(ANG2. LE.ANG1 )  ANG2=2 . 0*PI+ANG2 

FIND  ANGULAR  INCREMENT 

DI=  ABS ( FLOAT ( ISTP-IST  RT ) ) 

DJ=  ABS ( FLOAT (JSTP-JSTRT)) 

IINC=0 
JINC=0 

IF ( ISTRT . NE . ISTP )  IINC= ( 1STP-ISTRT ) /LABS ( ISTP-ISTRT ) 
IF( JSTRT . NE . JSTP)  JINC= (JSTP-JSTRT ) /TABS ( JSTP-JSTRT ) 
LAMDA=1 

IF( IINC . NE . 0 . AND . JINC . NE . 0 )  LAMDA=2 
DIFF=MAX1 (DI,DJ) 

ITER=DIFF-1 . 

IF ( LAMDA . EQ . 2 )  DIFF=2.*D1FF 
DELPHI=(ANG2-ANG1 )/DIFF 
WRITE (6,2008)  ANG 1 , ANG2 , DIFF , DELPHI 

CHECK  FOR  INPUT  ERROR 

IF ( LAMDA . NE . 2 . OR . DI . EQ . D J )  GO  TO  350 
WRITE (6, 200 3) 

GO  TO  150 
350  IO=ISTRT 
JO=JSTRT 
WRITE(6 ,2004) 

C 

C  INTERPOLATE 
C 

NPT=IABS(I2-I1 )+IABS( J2-J1 )-1 
DO  380  M=1 , ITER 

359  IF ( LAMDA . EQ . 2 )  GO  TO  360 
I=10+IINC 

J=JO+JINC 
CALL  MNIMX(I , J) 

NCODE(I ,  J)  =  1 

CALL  CIRCLE(ANG1 , DELPHI, RSTRT,ZSTRT,RC,ZC, I, J) 
WRxTE(6 ,2005)  I, J,AR(I, J) ,AZ(I , J) 

GO  TO  370 

360  I=IO+IINC 
J=JO 

NCODE(I , J)  =  1 

call  mnimx(i,j) 

CALL  CIRCLE ( ANG 1 , DELPHI , RSTRT , ZSTRT , RC , ZC , I , J ) 
WRITE(6 , 2005 )  I,J,AR(I,J) ,AZ(I, J) 

J=JO+JINC 
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NCO0E(1,J)=1 
CALL  MN1MX(1 , J) 

CALL  CIRCLE (ANG1 , DELPHI, RSTRT , ZSTRT , RC , ZC , I ,J) 

WR1TE(6,2005)  i,J,AR(I,J),AZ(I,J) 

370  10=1 
3S0  JO=J 

IF ( LAMDA .  NE .  2 )  CO  TO  390 
I=±0+IINC 
NCODE(I, J)=1 
call  mnimx(i,j) 

CALL  CIRCLE ( ANC 1 , DELPHI , RSTRT , ZSTRT , RC , ZC , I , J ) 

WRxTE(6,2005)  I,J,AR(I,J),AZ(I,J) 

390  IF(KAPPA.EQ.2)  GO  TO  150 
ISTRT=I2 
ISTP=I3 
JSTRT=J2 
JSTP=J3 
R3TRT=R2 
RSTPsflJ 
Z3TRT=Z2 
Z3TP=Z3 
KAPPA=2 

399  GO  TO  340 

C*  *********************************** 

C  CALCULATE  COORDINATES  OF  INTERIOR  POINTS 

C*  *********************************** 

400  IF(MAXJ.LE.2)  GO  TO  430 
J2=MAXJ-1 

DO  420  N= 1 ,500 
RESID=0 . 

DO  410  J=2,J2 
I 1 =IMIN( J)+1 
I2=IMAX( J)-1 
DO  410  1=11,12 

IF ( NCODE ( I , J ) . EQ . 1 )  GO  TO  410 

DR=(AR(I+1 , J)+AR(I-1 ,J)+AR(I,J+1)+AR(I,J-1))/4.-AR(I,J) 

DZ=(AZ(I+1 , J)+AZ(I-1 , J)+AZ( 1 , J+1 )+AZ (I , J— 1 ))/4. -AZ( I , J ) 
RESID=R£SID+AB3 ( DR ) +ABS ( DZ ) 

ARd,J)=ARd,J)  +  1  .  8*DR 
AZ(I,J)=AZ(I,J)+1 . 8*0Z 
410  CONTINUE 

IF(N.EQ.I)  RES1 =R£SID 
IF ( N . EQ . 1 . AND . RESID . EQ . 0 . ) GO  TO  430 
IF(RESID/RE31 .LT.1.E-5)  GO  TO  430 
420  CONTINUE 
430  WRITE(6 , 2009 )  N 

C*  *********************************** 

CALL  POINTS 

C*  *********************************** 

1000  FORMAT  (515) 
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1001  FORMAT  (3(213, 2F8. 3), 15) 

2000  FORMAT  (30H1  MESH  GENERATION  INFORMATION// 

1  41H0  MAXIMUM  VALUE  OF  I  IN  THE  MESH - 13/ 

2  41  HO  MAXIMUM  VALUE  OF  J  IN  THE  MESH - 13/ 

3  41H0  NUMBER  OF  LINE  SEGMENT  CARDS - 13/ 

4  41H0  NUMBER  OF  BOUNDARY  CONDITION  CARDS - 13/ 

5  41H0  NUMBER  OF  MATERIAL  BLOCK  CARDS - 13///) 

2001  FORMAT  (//88H  INPUT  II  J1  R1  Z1  12  J2  R2  Z 

12  13  J3  R3  Z3  IPTION/8X,3(2l4,2F8.4) ,16) 

2002  FORMAT  (5H  DI=F4.0,5H  DJ=F4.0,7H  DlFF=F4.0,7H  RINC=F8.3,7H  ZI 

1NC=F8.3,7H  ITER=I3,7H  IINC=I3,7H  JINC=I3,8H  KAPPA=I1) 

2003  FORMAT(1X,38H**BAD  INPUT— THIS  LINE  IS  NOT  DIAGONAL) 

2004  FORMAT  (30H  I  J  AR  AZ) 

2005  FORMAT  (215 , 2F1 1 .6) 

2006  FORMAT  (5 1 H  **  BAD  INPUT  -  THESE  POINTS  DO  NOT  DEFINE  A  CIRCLE,/, 
13X,6F12.4,10X,2E20.8) 

2007  F0RMAT(19H  CENTER  COORDINATE ,( F1 1 . 6 , 1 X , F1 1 . 6 , 1 X) ) 

2008  FORMAT  (7H  ANG1=F9.6,7H  ANG2=F9-6,7H  DIFF=F3.0,9H  DELPHI=F9.6) 

2009  FORMAT  (//30H  COORDINATES  CALCULATED  AFTER  I3,11H  ITERATIONS) 
RETURN 

END 

C 

C 

C  SUBROUTINE  MINV 

C 

C  PURPOSE 

C  INVERT  A  MATRIX 

C 

C  USAGE 

C  CALL  MINV(A,N,D,L,M) 

C 

C  DESCRIPTION  OF  PARAMETERS 

C  A  -  INPUT  MATRIX,  DESTROYED  IN  COMPUTATION  AND  REPLACED  BY 

C  RESULTANT  INVERSE. 

C  N  -  ORDER  OF  MATRIX  A 

C  D  -  RESULTANT  DETERMINANT 

C  L  -  WORK  VECTOR  OF  LENGTH  N 

C  M  -  WORK  VECTOR  OF  LENGTH  N 

C 

C  REMARKS 

C  MATRIX  A  MUST  BE  A  GENERAL  MATRIX 

C 

C  SUBROUTINES  AND  FUNCTION  SUBPROGRAMS  REQUIRED 

C  NONE 

C 

C  METHOD 

C  THE  STANDARD  GAUS3-J0RDAN  METHOD  IS  USED.  THE  DETERMINANT 

C  IS  ALSO  CALCULATED.  A  DETERMINANT  OF  ZERO  INDICATES  THAT 

C  THE  MATRIX  IS  SINGULAR. 

C 
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c 

SUBROUTINE  MiNV(A,N,D,L,M) 
dimension  a(i),l(i),m(i) 
c 

c  . 

c 

C  IF  A  DOUBLE  PKECiSION  VERSION  OF  THIS  ROUTINE  IS  DESIRbD,  THE 

C  C  IN  COLUMN  1  SHOULD  BE  REMOVED  FROM  THE  DOUBLE  PRECISION 

C  STATEMENT  WHICH  FOLLOWS. 

C 

C  DOUBLE  PRECISION  A,D,BIGA,HOLD 

C 

C  THE  C  MUST  ALSO  BE  REMOVED  FROM  DOUBLE  PRECISION  STATEMENTS 

C  APPEARING  lN  OTHER  ROUTINES  USED  IN  CONJUNCTION  WITH  THIS 

C  RoUTxNE. 

C 

C  THE  DOUBLE  PRECISION  VERSION  OF  THIS  SUBROUTINE  MUST  ALSO 

CONTAIN  DOUBLE  PRECISION  FORTRAN  FUNCTIONS.  ABS  IN  STATEMENT 
10  MUST  BE  CHANGED  TO  DABS. 


SEARCH  FOR  LARGEST  ELEMENT 

D=1.0 

NK=-N 

DO  80  K=  1 ,  N 

NK=NK+N 

L(K)=K 

M(K)=K 

KK=NK+K 

BIGA=A(KK) 

DO  20  J=K,N 
IZ=N*( J-1 ) 

DO  20  I=K , N 
IJ=IZ+I 

10  IF(ABS(BIGA)-ABS(A(IJ)))  15,20,20 
15  BiGA=A(IJ) 

L(K)=I 
M(K)=J 
20  CONTINUE 

INTERCHANGE  ROWS 

J=L(K) 

IF(J-K)  35,55,25 
25  KI=K-N 

DO  30  1=1, N 
KI=KI+N 
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HOLD=-A(KI ) 

JI=KI-K+J 
A(Ki)=A( JI) 

50  A( JI )  =HOLD 
C 

G  INTERCHANGE  COLUMNS 

C 

55  I=M(K) 

IF(I-K)  45,45, 38 
38  JP=N*( 1-1 ) 

DO  40  J= 1 , N 
JK=NK+J 
Jl=JP+J 
H0LD=-A( JK) 

A( JK)=A( JI) 

40  A( JI)  =H0LD 

DIVIDE  COLUMN  BI  MINUS  PIVOT  (VALUE  UF  PIVOT  ELEMENT  IS 
CONTAINED  IN  BIGA) 

45  IF (BIGA)  43,46,48 

46  D=0 . 0 
RETURN 

48  DO  55  1=1, N 

IFU-K)  50,55,50 
50  IK=NK+I 

A( IK)=A( IK)/( -BIGA) 

55  CONTINUE 

REDUCE  MATRIX 

DO  65  1=1, N 
IK=NK+i 
H0LD=A( IK) 

IJ=I-N 
DO  65  J=1,N 
IJ=I J+N 

IF(I-K)  60,65,60 
60  IF(J-K)  62,65,62 
62  KJ=IJ-I+K 

A(IJ)=HOLD*A(KJ)+A(IJ) 

65  CONTINUE 

DIVIDE  ROW  BY  PIVOT 

KJ=K-N 
DO  75  J=1,N 
KJ=KJ+N 

IF(J-K)  70,75,70 
70  A ( K J ) = A ( K J ) /BIGA 
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75  CONTINUE 

product  of  pivots 


D=D*BIGA 


C  REPLACE  PIVOT  BY  RECIPROCAL 

C 

A(KK)=1.0/B±GA 
80  CONTINUE 

C 

C  FINAL  ROW  AND  COLUMN  INTERCHANGE 

C 

K=N 

100  K= (K-1 ) 

IF(K)  150,150,105 
105  I=L(K) 

IF(i-K)  120,120,103 
108  JQ=N*(K-1 ) 

Jfi=N*( 1-1  ) 

DO  110  J=1 , N 
JK=JQ+J 
HOLD=A( JK) 

JI=JR+J 
A( JK)=-A( JI) 

110  A( JI)  =HOLD 
120  J=M(K) 

IF(J-K)  100,100,125 
125  K1=K-N 

DO  130  1=1, N 
KI=KI+N 
HOLD=A(KI ) 

JI=KI-iC+J 
A(K1)=-A( JI) 

130  A(J1)  =HOLD 
GO  TO  100 
150  RETURN 
END 

SUBROUTINE  MNIMX(I,J) 

COMMON/TD/ IMIN (100), IMAX( 100), JMIN(25 ) , JMAX(25 ) , MAXI ,MAXJ ,NMTL , NBC 

IF( J . LT . JMIN (I) )  JMIN(I)=J 

IF( J.GT . JMAX(I) )  JMAX(I)=J 

IF ( I . LT . IMIN ( J ) )  IMIN(J)=I 

1F(I.GT.IMAX(J))  IMAX( J)=I 

RETURN 

END 

SUBROUTINE  MODIFY (NEQ, N , U) 

COMMON/SOLVE/B( 1 62 ) , A( 1 62 , 81 ) , NUMBLK , MBAND 

DO  10  M=2, MBAND 

K=N-M+1 
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IF(K.LE.O)  GO  TO  5 
B(K)=B(K)-A(K,M)*U 
A(K,M)=0.00 
5  K=N+M-1 

lF(NEQ.LT.K)  GO  TO  10 
B(K)=B(K)-A(N ,M) *0 
A(N ,M)=0.00 
10  CONTINUE 
A(N,1)=1 .00 
B(N)=U 
RETURN 
END 

SUBROUTINE  MPLOT 
INTEGER  CODE 

COMMON/TD/ IMIN (100), IMaX( 100 ) , JMIN (25 ) , JMAX(25 ) ,MAXI ,MAXJ ,NMTL , NBC 
COMMON /NPOATA/fl( 1000) , CODEC  1000) ,XR( 1000) ,Z( 1000) ,XZ(1000) , 
1NPNUMC25, 80 ) ,T( 1000) ,XT( 1000) 

REAL  X(100) ,Y( 100) ,TX(2) ,TY(2) ,TITLE(20) ,ZMAX 
READ  (5,1000)  TITLE , RMAX , ZMAX 
C  CALL  CCP2SY  (0.7, 0.2, 0.2, TITLE, 0.0,80) 

C  CALL  CCP1PL  (0.7, 0.7, -3) 

TX( 1 )=0.00 
TY(1 )=0.00 
TX(2)=RMAX/9 . 0 
TY(2)=RMAX/9.0 
ZMAX=ZMAX*TY (2 )+2 . 0 
IF  (ZMAX. LT. 17.0)  ZMAX=17.0 
DO  100  J=1 ,MAXJ 
NSTART=IMIN( J) 

N3T0P=IMAX( J) 

N  =  0 

DO  101  I=NSTART , NSTOP 
N=N+1 

NP=NPNUM(I, J) 

Y(N  )  =  R(NP ) 

101  X(N)=Z(NP) 

C  CALL  CCP6LN  (X,Y,N, 1 ,TX,TY) 

100  CONTINUE 

DO  102  1=1, MAXI 
NSTART=JMIN(I) 

N3T0P=JMAX(I) 

N=0 

DO  103  J=NSTART, NSTOP 
N=N+1 

NP=NPNUM(I, J) 

Y(N)=R(NP) 

103  X(N)=Z(NP) 

C  CALL  CCP6LN  (X,Y,N, 1 ,TX,TY) 

102  CONTINUE 

C  CALL  CCP1PL  (ZMAX, -0.7, -3) 

101 


1000  FORMAT  (20A4/2F 10 . 0 ) 

RETURN 

END 

SUBROUTINE  NAXSTF( II , JJ , KK) 

INTEGER  CODE 

COMMON /MATP/RO (6 ) ,E(12,l6,6),EE(l6) , AOFTS(6) 

COMMON  /  BASIC/ ACELZ ,  ANGVEL ,  ANGACC ,  TREF ,  VOL ,  NUMNP ,  NUMEL ,  NUMPC ,  NUMSC , 
1NUMST 

COMMON/ARG/RRR(5),ZZZ(5),RR(4),ZZ(4),S(15,15),P(15),TT(6), 

1 H ( 6 , 15) ,CRZ(6 , 6) ,XI( 10 ) , ANGLE(4) ,Slu( 18) , EPS( 1 8) ,N 
COMMON/NPDATA/R( 1000 ) , CODE ( 1 000 ) ,XR( 1000) ,Z( 1000 ), XZ( 1000) , 
1NPNUM(25, 80), T( 1000), XT( 1000) 

COMMON/ELDATA/BETA(1000) , EPR( 1 000 ) , PR(200 ) , SH (200 ) , IX( 1 000 , 5 ) , 

IIP (200) , JP(200) ,IS(200) ,  JS(200 ) , ALPHA ( 1000) , IT (200) , JT(200 ) , 
2ST(200 ) 

C0MM0N/NXQUAD/AR1 

C01vlM0N/N0NAXI/S  1  ( 30 , 30 ) ,  P 1  ( 30 ) ,  THETA , BS 1 ( 6 , 30 ) 

DIMENSION  C(18,18),B(18, 18),B1(6,18),B2(6,18),B3(6,18),B4(6,18), 

1  B5(6,18),B6(6,18),B1A(6,18),B1B(6,18) ,B2A(6, 18) ,B2B(6, 18 

2  ) ,B3A(6,l8),B3B(6,l8),B4A(6,l8),B4B(6,l8),B5A(6, 18), 

3  B5B(6, 18) ,B6A(6, 18) ,B6B(6, 18) 

C  ZERO  MATRICES 

DO  100  1=1,18 

DO  100  J= 1  , 1 8 
100  C(I,J)=  0.0 
DO  110  1=1,6 
DO  110  J=  1  , 1 8 
B1 (I , J)  =0.0 
B2(I,J)  =0.0 
B3(I,J)  =0.0 
B4(I,J)  =0.0 
B5(i,J)  =0.0 
110  B6(I,J)  =0.0 
RR( 1 )  =  RRR(II) 

RR(2)  =  RRR(JJ) 

RR(3)  =  RRR(KK) 

ZZ(1)  =  ZZZ(II) 

ZZ(2)  =  ZZZ(JJ) 

ZZ(3)  =  ZZZ(KK) 

C0MM=RR(2 )*( ZZ(3 )-ZZ( 1 ) )+RR ( 1 )*(ZZ(2)-ZZ( 3) )+RR(3 )*( ZZ( 1 )-ZZ(2 ) ) 

C  FILL  C  INVERSE 

C(  1  , 1  )=  (  RR(2)*ZZ(3)  -RR(3)*  ZZ(2))  /  COMM 

C(  1 , 4)=  (  RR(3 )*ZZ(  1 )  -RR(1 )*  ZZ(3) )  /  COMM 

C( 1 ,7)=  (  Rfl(1)«ZZ(2)  -RR(2)*  ZZ ( 1 ) )  /  COMM 

C(2 , 1 ) =  (  ZZ(2)  -  ZZ( 3 ) )  /  COMM 

C(2 , 4)=  (  ZZ(3)  -  ZZ(1))  /  COMM 

C(2,7)=  (  ZZ(1)  -  ZZ(2) )  /  COMM 

C(3 , 1 )=  (  RR(3)  -  RR(2) )  /  COMM 

0(3,4)=  (  RR ( 1 )  -  RR(3) )  /  COMM 

C(3,7)=  (  RR(2)  -  RR ( 1 ) )  /  COMM 
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C(4,2)=  C( 1 , 1 ) 

C(4,5)=  C( 1 ,4) 

C ( 4 , 8 )  =  C(1, 7) 

C(5,2)=  C(2 , 1 ) 

C(5,5)=  0(2,4) 

C(5,8)=  C(2 , 7 ) 

C(6,2)=  C( 8 , 1 ) 

C(6 , 5)=  C(3, 4) 

C ( 6 , 8 ) =  C(S , 7 ) 

C ( 7 , 3 ) =  C(1,1) 

C ( 7 , 6 ) =  C(1, 4) 

C(7,9)=  C(1, 7) 

C(8 , 3)=  C(2 , 1 ) 

C(8 , 6)=  C(2 , 4) 

C(8,9)=  C ( 2 , 7 ) 

C(9,3)=  C(3, 1 ) 

C(9,6)=  C(3,4) 

C(9 ,9)  =  C(3,7) 

DO  120  1=10,18 
DO  120  J=  1,9 
II  =  1-9 
J1=J+9 

C(I,J)  =(-1. /THETA)  *  C ( I 1 , J ) 
C(I,J1)=(  1. /THETA)  *  C(I1,J) 
120  CONTINUE 
C  FILL  B  MATRICES 
C  B1  CONSTANT  TERMS 
C  B2  THETA  TERMS 
C  B3  1/R  TERMS 
C  B4  THETA/  R  TERMS 
C  B5  Z/R  TERMS 
C  B6  THETA  *Z/R  TERMS 
DO  130  J=  1,1 8 
B 1  ( 1 , J )  =  C(2 , J) 

B1 (2 , J)  =  C(6, J) 

B1 (3, J)  =  C(2,J)+C(17,J) 
B1(4,J)  =  C(3,J)+C(5,J) 

B1 (5, J)  =  C(9 , J)  +C ( 1 4 , J ) 

B1 (6 , J)  =  C( 1 1 , J) 

B2 ( 1 , J )  =  C( 1 1 , J) 

B2(2, J)  =  C( 15, J) 

B2(3, J)  =  C( 1 1 , J) 

B2(4, J)  =  C(12,J)+C(14,J) 

B2(5 , J)  =  C ( 1 8 , J ) 

B3(3,J)  =  C(1,J)+  C ( 1 6 , J ) 
B3(5,J)  =  C( 13, J) 

B3(6,J)  =  C( 10 , J)  -  C(7 , J) 
B4(3,J)  =  C(10,J) 

B4(6 , J)  =  -C( 16, J) 

B5(3, J)  =  C(3,J)  +C ( 1 8 , J ) 
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C( 15 , J) 

C(  12 , J)-C(9 , J) 
C( 12 , J) 

— C  ( 1 8 ,  J) 

BT  *  D  *  B 


B5(5,J)  = 

B5(6,J)  = 

66(3, J)  = 

B6(6 , J)  = 

150  CONTINUE 
NOW  CALCULATE 
CALL  INTER 

THETA2  =  (THETA  **2)/2.0 
THETA3  =  (THETA  **3)/3.0 
DO  140  1=1,6 
DO  140  J=1 ,18 
B1A(i,J)=(B1 (I, J)*XI(1 ) 

1  (B2(I,J)*XI(1) 

B2A (I,J)=(B1 (I, J) *XI ( 1 ) 

1  +  (B2 (I , J) *XI ( 1 ) 

B3A(I,J)=(B1(I,J)*XI(2) 

1  +(B2(I,J)*XI(2) 

B4A(I,J)=(B1(I,J)«XI(2) 

1  +(B2(I,J)«XI(2) 

B5A(I,J)=(B1(I,J)*XI(4) 

1  +  (B2(I, J)*XI(4) 

B6A(I,J)=(B1(I,J)«XI(4) 

1  +  (B2(I,J)*XI(4) 

140  CONTINUE 

DO  150  1=1,6 
DO-  150  K=  1 ,  18 
B1B(I , K) =  0.0 


+B3 ( I , J ) * 

XI  (2) 

+ 

B5(I, J)» 

XI (4 ) )  * 

THETA 

+B4(I, J)* 

XI  (2 ) 

+ 

B6(I, J)* 

XI (4 ) ) * 

THETA2 

+B3(I,J)« 

XI  (2 ) 

+ 

B5(i,J)» 

XI (4 ) ) * 

THETA2 

+B4(I, J)* 

XI  (2) 

+ 

B6(I, J)* 

XI  ( 4 ) ) * 

THETAj 

+B3 ( I , J ) * 

XI  (3) 

+ 

B5(I,J)» 

Xx(5))* 

theta 

+B4(I,J)» 

XI(3) 

+ 

B6 (I , J) * 

XI(5) )* 

THETA2 

+B3(1,J)* 

XI(3) 

+ 

B5(I , J) * 

XI(5) )* 

THETA2 

+B4(I,J)« 

xi  (3) 

+ 

B6 (I , J) * 

XI (5 ) ) * 

THETA 3 

+B3(I,J)« 

xi(5) 

+ 

B5 (1 , J )  * 

XI(6))» 

THETA 

+B4(I, J)« 

XI  (5 ) 

+ 

B6(I,J)« 

XI ( 6 ) ) * 

THETA2 

+B3(I,J)« 

XI(5) 

+ 

B5(I,J)« 

XI(6) )* 

THETA2 

+B4(I, J)« 

XI  (5 ) 

+ 

B6(I,J)« 

XI(6 ) ) * 

THETA 3 

B2B(I,K)= 
B3B(I,K)= 
B4B(I,K)= 
B5B(I,K)= 
B6B(I,K)= 
DO  150 


0.0 
0.0 
0.0 
0.0 
0.0 
J=1 , 6 

B1B(I,K) 
B2B(x,K) 
B3B(I,K) 
B4B(I,K) 
B5B(I,K) 
B6B(I,K) 


B1A(J,K) 
B2A( J , K) 
B3A(J,K) 
B4A(J,K) 
B5A(J,K) 
B6A(J,K) 


B1B(I,K)  =  B1B(I,K)  +  CRZ(I,J) 

B2B( I , K)  =  B2B(x,K)  +  CRZ(I,J) 

B3B(I,K)  =  B3B(I,K)  +  CRZ(I, J) 

B4B( I , K)  =  B4B(I,K)  +  CRZ( I , J ) 

B5B(I,K)  =  B5B(i,K)  +  CRZ(I,J) 

B6B( I , K)  =  B6B( I , K)  +  CRZ(I,J) 

1^0  CONTINUE 

DO  160  1=1,18 

DO  160  K= 1 , 1 8 
B(I,K)=0 .0 
DO  160  J  =  1 , 6 

B(I,K)  =  B(I,K)  +  B1(J,I)»  B1B( J , K)+B2( J , I) *B2B( J , K)+B3 ( J , I)* 

1  B3B(J,K)+B4(J,I)*B4B(J,K)+B5(J,I)*B5B(J,K)+B6(J,I)*B6B(J,K) 

1 60  CONTINUE 
250  CONTINUE 

C  B(I,K)  NOW  CONTAINS  THE  STIFFNESS  MATRIX  FOR  ONE  TRIANGULAR  ELEMENT 
AR1  =  AR1  +  XI( 1 )  *THETA 
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255 


DO  235  K=1 ,6 
DO  235  1=1,3 
8S1(K,3*iI-3+D  = 
BS1(K,3*JJ-3+D  = 
B31(K,3*KK-3+I)  = 
BS1(K, 3*11+1+12)= 
B31(K,3*JJ+I+12)= 
831 (K , 3*KK+I+12)= 
CONTINUE 


B31(K, 3*11-3+1)  +B 1 A ( K , I  ) 
BS1(K,3*JJ-3+I)  +B1A(K , l+j  ) 
B31 (K , 3*KK-3+I)  +B1A(K , 1+6  ) 
B31(K,3*H+12+I)+B1A(K,I+9  ') 
B31 (K, 3*JJ+1 2+1 )+B1A(K, 1+12) 
BS1 (K, 3*KK+12+I)+B1A(K, 1+15) 


182 


180 


1IM  =  3*  Ii  -3 

JJM  =  3*  JJ  -3 

KKM  =  3*  KK  -3 

DO  170  K= 1 

,4 

DO  170  1=1 

,3 

DO  170  J=1 

,3 

IF(K.EQ. 1 

.OR. 

K.EQ.2) 

11=1 

IF(K. EQ. 3 

.OR. 

K.EQ.4) 

11=1  +9 

1F(K.EQ. 1 

.OR. 

K.EQ.3) 

J1=J 

IF(K.EQ.2 

.OR. 

K.EQ.4) 

J1=J 

IF(K.EQ.1 

.OR. 

K.EQ.2) 

K1  =0 

IF(K.EQ.3 

.OR. 

K.EQ.4) 

K 1  =  1 5 

iF(K.EQ. 1 

.OR. 

K.EQ.3) 

PS 

IV 

II 

o 

IF(K . EQ. 2 

.OR. 

K.EQ.4) 

t_n 

n 

cv 

+9 


170 


KK2=KKM 
II2=1IM 
JJ2=JJM 
KK 1 =KKM 
JJ1=JJM 
II1=IIM 

SKII1+I+K1  ,Ii2+J+K2) 

31 (1I1+I+K1 , JJ2+J+K2) 

31 (1I1+I+K1 , KK2+J+K2) 

31 ( JJ 1+I+K1 , II2+J+K2) 

31 (JJ1+I+K1 , JJ2+J+K2) 

3 1 ( J J 1 +I+K 1 , KK2+J+K2 ) 

31 (KK1+1+K1 , II2+J+K2) 

SI (KK1+I+K1 , JJ2+J+K2) 

31 (KK1+I+K1 ,KK2+J+K2) 

CONTINUE 

RETURN 

END 

FUNCTION  NODE(l , J) 

COMMON /TD/IMIN (100), IMAX( 100), JMIN ( 25 ) , JMAX( 25 ) ,MAXI , MAXJ , NMTL , NBC 
N0DE=0 

DO  100  JJ=1  ,J 
NSTART=IMIN ( JJ ) 

NSTOP=IMAX( JJ ) 

DO  100  II=N3TART , NSTOP 
NODE=NODE+1 

IF( JJ . EQ. J . AND . II .EQ.I)  RETURN 


SK1I1+I+K1  ,II2+J+K2) 
SKII1+I+K1 ,  JJ2+J+K2) 
S1(II1+I+K1 , KK2+J+K2) 
SI (JJ1+I+K1 ,II2+J+K2) 
SKJJ1+I+K1  ,JJ2+J+K2) 
SKJJ1+I+K1 ,  KK2+J+K2) 
SI (KK1+I+K1 ,II2+J+K2) 
SKKK1+I+K1  ,JJ2+J+K2) 
31 (KK1+1+K1 , KK2+J+K2) 


+B( II , J 1 ) 

+B(I1 , J1+3) 
+B(I1 , J1+6) 
+B(I1+3, J1 ) 
+B(I1+3, J1+3) 
+B(I1+3, J1+6) 

+  B(I1+6,J1) 

+  B(I1+6, J1+3) 
+  B(I1+6, J1+6) 
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100  CONTINUE 
RETURN 
END 

SUBROUTINE  POINTS 
INTEGER  CODE 

COMMON /BASIC/ ACELZ , ANG VEL , ANGACC , TREF , VOL , NUMNP , NJMEL , NUMPC , NUMSC , 
1N0M3T 

COMMON/MATP/RO(6 ) ,E(12, 1 6 , 6 ) ,E£( 16 ) , AOFTS(6) 

COMMON /NPDATA/ R ( 1 000 ) , CODE (1000) , XR( 1 000 ) , Z(  1 000 ) , XZ( 1 000 ) , 
1NPNUM(25,80) ,T( 1000) ,Xf( 1000) 

COMMON/ELDATA/BETA(  1000)  ,EPR(1000)  ,PR(200)  ,SH(200) , IX(1000,5) , 
1IP(200)  ,JP(200) , IS (200) ,JS(200) , ALPHA ( 1000) , IT (200) ,JT(200) , 
2ST(200 ) 

COMMON/ SOLVE/X(4428 ) , Y ( 4428 ) , TEM( 4428 ) , NdMTC ,MBAND 

COMMON/TD/ IMIN (100), IMAX( 100), JMiN (25 ) , JMAX(25 ) ,MaX1 ,MAXJ , NMTL , NBC 

COMMON/PLANE/NPP 

DIMENSION  AR ( 25 , 80 ) , AZ ( 25 , 80 ) , MATRIL( 100,5), BLKANG (100), BLKALF ( 1 
100) 

DIMENSION  IBNG(IOO) ,NBNG(100) 

EQUIVALENCE  (R( 1 ) ,AR) , (Z( 1 ) , AZ) 

C  ESTABLISH  NODAL  POINT  INFORMATION 

<Z*  *********************************** 

NEL=0 
N0D3UM=0 
DO  100  J= 1 ,MAXJ 
NSTART=IMIN(J) 

NSTOP=IMAX( J) 

DO  100  ±=N3TART , NSTOP 
100  N0DSUM=N0DSUM+1 
NEL3UM=0 
JJMAX=MAXJ-1 
DO  110  JJ=1 , JJMAX 
N3TOP=MINO(IMAX(JJ) , IMAX( JJ+1 ) ) -1 
NST ART =MAXO (IMIN (JJ ) , IMiN ( JJ+1 ) ) 

DO  110  II=NSTART, NSTOP 
110  NELSUM=NEL3UM+1 
NUMNP=N0D3UM 
NUMEL=NELSUM 
DO  120  J= 1 ,MAXJ 
N3TART=IMIN(J) 

NSTOP=IMAX( J) 

DO  120  I=NSTART , NSTOP 
NPNUM(I, J)=NODE(I, J) 

NP=NPNUM( I , J) 

R(NP)=AR(I , J) 

120  Z(NP)=AZ(I, J) 

C*  *********************************** 

C  READ  AND  ASSIGN  BOUNDARY  CONDITIONS 

C*  *********************************** 

C  INITIALIZE 
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^3*  ft********************************** 

DO  130  1=1 ,NUMNP 
CODE(I)=0 

IF(R(I) .EQ.O. .AND.NPP.EQ.O)  C0UE(I)=1 . 

XR(I)=0 . 

XZ(1)=0. 

XT ( I )  =  0 . 0 
130  T(1)=0. 

IF (NBC . EQ. 0 )  GO  TO  210 
DO  200  IBCON=1 , NBC 

READ(5, 1002)  II ,12, J1 , J2 , ICN , RCON , ZCON , TCON 
DO  200  1=11,12 
DO  200  J=J1,J2 
NP=NPNUM(I,J) 

CODE (NP)= ICN 
XR(NP)=RCON 
XT(NP)=TCON 
200  XZ(NP)=ZCON 
210  MPRINT=0 

DO  230  J=1 ,MAXJ 
NSTART=IMIN ( J ) 

N3T0P=IMAX( J) 

DO  230  I=N3TART ,NST0P 
NP=NPNUM(I, J) 

IF ( MPRINT . NE . 0 )  GO  TO  220 
WRITE(6 ,2000) 

MPRINT=59 

220  MPRINT=MPRINT-1 

230  WRITE(6,2001 )  I, J,NP,CODE(NP) ,R(NP) ,Z(NP) ,XR(NP) ,XZ(NP) ,XT(NP) 

(3#####*#**#******#########*##**#***** 

C  ASSIGN  MATERIALS  IN  BLOCKS 

C*  *********************************** 

DO  300  Ml =1 ,NJMEL 
300  ±X(M1 ,5)=0 

DO  310  IMTL=1 ,NMTL 

READ  (5,1000)  MTL , ( MATRIL ( IMTL , IM) , IM=2 , 5 ) , BLKANG ( IMTL ) , BLKALF ( IMT 
1 L ) , IBNG ( IMTL ) , NBNG ( IMTL ) 

310  MATRIL (IMTL, 1)=MTL 

Q*  *********************************** 

C  ESTABLISH  ELEMENT  INFORMATION 

C*  *********************************** 

JJMAX=MAXJ-1 
N=0 
MTL  =  1 
KTL=  1 

DO  440  JJ=1 , JJMAX 
NSTOP=MINO(IMAX( JJ ) , IMAX( JJ+1 ) )-1 
NSTAfiT=MAXO ( IMIN ( J J ) , IMIN ( J J+1 ) ) 

DO  440  II=NSTART , N3T0P 
NEL=NEL+1 
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DO  400  1MTL= 1 , NMTL 
IF(Ii.LT .MATR1L(IMTL , 2) )  GO  TO  400 
iFdI.GE.MATHIL(IMTL,3))  GO  TO  400 
IF(JJ.LT .MATfiIL(IMTLf 4))  GO  TO  400 
IF( JJ . GE .MATRIL(IMTL , 5) )  GO  TO  400 
KAT=IMTL 

MaT=MATRIL(IMTL,1) 

400  CONTINUE 

IF(KAT.EQ.KTL)  GO  TO  410 

KTL=KAT 

MTL=MAT 

GO  TO  420 

410  IF(II.EQ.NSTART)  GO  TO  420 

IF ( JJ.NE.JJMAX. OR. II. NE.NSTOP)  GO  TO  440 

M=NEL+1 

IANG=ICNG 

NANG=NCNG 

GO  TO  421 

420  l=NPNUMdI, JJ) 

J=I+1 

K=NPNUM(II+1 ,JJ+1 ) 

L=K-1 

M=NEL 

IX(M, 1 )=I 

1X(M,2)=J 

IX(M,3)=K 

IX(M, 4)=L 

IX(M, 5)=MTL 

BETA ( M) =BLKANG (KTL ) 

ALPHA (M)=BLKALF(KTL) 

IANG=ICNG 

NANG=NCNG 

ICNG=IBNG(KTL) 

NCNG=NBNG(KTL) 

421  NC=2 
430  N=N+1 

IF(M.LE.N)  GO  TO  440 
IX(N , 1 )=IX(N-1 , 1 )+1 
IX(N,2)=IX(N-1 ,2)+1 
IX(N,3)=IX(N-1 ,3)+1 
IX(N,4)=IX(N-1 ,4)+1 
IX(N,5)=IX(N-1 ,5) 

8ETA(N)=BETA(N-1 ) 

IF(IANG.EQ.I)  GO  TO  442 
ALPHA(N )=ALPHA(N-1 ) 

GO  TO  443 
442  CONTINUE 

IF(NC.GT.NANG)  GO  TO  444 
ALPHA(N)=ALPHA(N-1 ) 

GO  TO  443 
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o  o 


444  NC=1 

ALPHA (N)=-ALPHA(N-1 ) 

443  CONTINUE 
NC=NC+1 

IF(M.GT.N)  GO  TO  430 
440  CONTINUE 

IF(NUMNP.GT .2000)  WRiTE (6, 2002) 

it*********************************** 

SET  NODAL  POINT  TEMPERATURE  TO  REFERENCE  TEMPERATURE 

^tt*********************************** 

IF(NUMTC.NE.O)  RETURN 
DO  500  N=1 ,NUMNP 
500  T (N )=TR£F 

1000  FORMAT  (515, 2F10. 0,215) 

1002  FORMAT (415 , 110, 3F 10.0) 

2000  FORMAT  (128H1  I  J  NP  TYPE  R-ORDINATE  Z-ORDINA 

1TE  R  LOAD  OR  DISPLACEMENT  Z  LOAD  OR  DISPLACEMENT  T  LOAD  OR  DISP 
2LACEMENT) 

2001  FORMAT  (215, 16 , 112 ,F13 . 6 ,F14 . 6 ,£26 . 7 , E24 . 7 ,E24 . 7) 

2002  FORMAT  (35H  BAD  INPUT  -  TOO  MANY  NODAL  POINTS) 

RETURN 

END 

SUBROUTINE  QUAD 
INTEGER  CODE 

REAL  NU3N , NUTN , NUTS , NUNS , NUNT , NUST 
DIMENSION  DUMMY (6, 6) , DUMMY1 (6,6) 

COMMON /BAS1C/ACELZ , ANGVEL , ANGACC , TREF , VOL , NUMNP , NUMEL , NUMPC , NUMSC , 
1NUMST 

COMMON /N  XQUAD/AR 1 

COMMON /N0NAX1/S 1 ( 30 , 30 ) , P 1 ( 30 ) , THETA , BS 1 ( 6 , 30 ) 

COMMON /MATP/R0(6) ,E(12,l6,6) ,EE(16) ,A0FTS(6) 

COMMON /NPDATA/R( 1000) , CODE ( 1000 ) ,XR ( 1 000 ), Z( 1 000 ) ,XZ(1000) , 
1NPNUM(25,80) ,T( 1000) , XT (1000) 

COMMON /ELDATA/BETA( 1000 ) ,EPR( 1000 ) , PR(200) ,SH(200) ,IX( 1000,5) , 

IIP (200)  ,JP(200)  ,13(200) , JS(200) , ALPHA ( 1000) , IT (200) , JT(200) , 
2ST(200 ) 

C0MM0N/ARG/RRR(5) ,ZZZ(5),RR(4),ZZ(4) ,S(15, 15) ,P(15),TT(6) , 

1 H ( 6 , 1 5 )  ,  CRZ ( 6 , 6 ) , XI ( 1 0 ) , ANGLE ( 4 ) , SIG ( 1 8 ) , EPS ( 1 8 ) , N 
C0MM0N/RESULT/B3(6 , 15),D(6,6),C(6,6) ,AR,BB(6,9) ,CNS(6,6) 

COMMON /PLANE/NPP 

COMMON/DUM1/S1TEM(3,30) ,S1T(24,24) ,TS(6,24) 

DIMENSION  S2T(24,6) 

DIMENSION  BS1T(6 , 3)  ,P1T(3) 

II =IX(N , 1 ) 

J1=IX(N,2) 

K1 =IX(N ,3) 

L1=IX(N ,4) 

MTYPE=IX(N ,5) 

IX(N , 5 )=-IX(N ,5) 

C*  *********************************** 


109 


c  interpolate;  material  properties 

C*  *********************************** 

DO  100  1=1  ,  12 
100  EE(I)=E(1 ,1+1 , MTIPE) 

DO  110  1=1,6 
DO  110  J=1 ,6 
CNS ( I , J ) =0 . 00 
C(I,J)=0.00 
110  D(±,J)=0.00 

C*  if********************************** 

C  FORM  STRESS-STRAIN  RELATIONSHIP  IN  N-S-T  SYSTEM 

C*  *********************************** 

NUNS=EE(4 ) 

NUNT=E£(5) 

NUST=EE(6 ) 

NUSN= (E£(2 )*NUN3)/EE( 1 ) 

NUTN=(EE(i)*NUNT)/EE(1 ) 

NUTS=  (EE  ( 3 )  *N(JST )  /EE  (2 ) 

D1V= 1 . 00-NUN3*NUSN-NUST*NUTS-NUNT*NUTN-NUSN*NUNT*NUT3 
1 -NUN5*NUTN*NUST 

CNS(1 , 1 )=EE( 1 )*( 1 . 00-NUST*N UTS) /DIV 
CNS ( 1 , 2 ) =E£ ( 2 ) * ( NUNS+NUNT*NUTS ) /DI V 
CNS ( 1 , 3 ) =££ (3 ) * (NUNT+NUN3*NUST ) /DIV 
CNS(2 , 1 )=CN3( 1,2) 

CNS ( 2 , 2 ) =EE ( 2 ) * ( 1 . 00-NUNT *NUTN ) /DIV 
CNS ( 2 , 3 ) =EE ( 3 ) * ( NUST+NUSN  *NUNT ) /DIV 
CNS (3 , 1 )=CNS( 1,3) 

CN3(3,2)=CNS(2,3) 

CNS( j ,3)=EE(3)*(1. 00-NUNS*NUSN ) /DIV 
CNS (4 ,4)=EE(7 ) 

CNS(5,5)=EE(8) 

CNS(6 , 6)=£E(9 ) 

C  SET  UP  STRAIN  TRANSFORM  TO  N-S-T  SYSTEM 
SINA=SIN (ALPHA (N)) 

COSA=COS(ALPHA(N) ) 

S2=SINA**2 
C2=C0SA**2 
SC=SINA*COSA 
D( 1 , 1 ) =C2 
D ( 1 , 3 ) =S2 
D(1 ,6)=-SC 
D(2 , 1 )=S2 
D(2 , 3)  =  C2 
D(2 , 6)=SC 
D(3,2)=1 .00 
D(4, 1 )=2.00*SC 
D(4,3)=-2.00*SC 
D(4,6)=C2-S2 
D(5,4)=SINA 
D(5 , 5)=COSA 


no 


D(6 , 4)=C0SA 
D(6,5)=-SINA 

C  SET  UP  STRAIN  TRANSFORMATION  TO  R-Z-T  SYSTEM 
SINB=SIN (BETA(N ) ) 

COSB=COS(BETA(N) ) 

S2=SINB**2 
C2=COSB**2 
SC=SINB*COSB 
C( 1 , 1 )=S2 
C(1 ,2)=C2 
C(1 ,4)=SC 
C(2 , 1 )=C2 
C ( 2 ,2)=S2 
C(2,4)=-SC 
C(3,3)=1.00 
C ( 4 , 1 )=-2 . 00*SC 
C ( 4 ,2)=2.00*SC 
C ( 4 , 4) =S2-C2 
C(5,5)=SINB 
0(5 , 6)=-C0SB 
C(6,5)=COSB 
C(6,6)=SINB 

C  CALCULATE  CRZ  MATRIX 
DO  120  1=1,6 
DO  120  J  =  1 , 6 
DUMMYU,  J)=0.00 
DO  120  K=1 ,6 

120  DUMMY ( I , J ) =DUMMY (I,J)+D(I,K)*C(K,J) 

DO  130  1=1,6 
DO  130  J=1 , 6 
DUMMY 1 (I, J)=0 .00 
DO  130  K=1 ,6 

130  DUMMY 1 (I, J)=DUMMY1 (I, J)+CNS(I,K)*DUMMY(K, J) 

DO  140  1=1,6 
DO  140  J=1 ,6 
DUMMY  ( I ,  J )  =  0 . 00 
DO  140  K=1 , 6 

140  DUMMY ( I , J ) =DUMMY (I,J)+D(K,I) *DUMMY 1 (K , J ) 

DO  160  1=1,6 
DO  150  J  =  1 , 6 
CRZ(I, J)=0.00 
DO  150  K=1 , 6 

150  CRZ(I, J)=CRZ(I, J)+C(K,I)*DUMMY(K, J) 

TT(I)=0.00 
DO  160  M=1 ,6 
P (M) =0 . 00 
DO  161  11=1,3 

IF(AOFTS(MTYPE) .EQ. 1 .)  P(M)=CNS(M,II)»EE(II+9) 

161  P(M)=P(M)+(T(N)-TREF)*CNS(M,II)*EE(II+9) 

DO  160  K=1 , 6 
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160  TT(I)=TT(I)+C(K,I)*D(M,K)»P(M) 


FORM  QUADRILATERAL  STIFFNESS  MATRIX 
RRR(5)=(R(I1)+R(J1)+R(K1 )+R(L1))/4. 
ZZZ(5)=(Z(I1)+Z(J1 )+Z(K1 )+Z(L1))/4. 

DO  200  M=1,4 
MM=1X(N,M) 

IF(NPP.NE.O)  GO  TO  190 

IF(R(MM) . EQ. 0 . .AND. CODE (MM) .EQ.O. ) CODE (MM) 
190  RRR(M)=R(MM) 

200  ZZZ(M)=Z(MM) 

DO  220  11=1 ,15 
P1(II)=0.0 
PI (11+15)  =0.0 
P(II)=0.00 
DO  220  JJ=1 , 15 
220  S(II,JJ)=0.00 
V0L=0 . 

DO  90  1=1,6 
DO  90  J= 1 ,15 
BS1(I,J)=0.0 
BS1 (I, J+15)  =  0.0 

90  BS(I, J)=0.00 
AR=0 . 00 

240  CALL  TRISTF(4, 1 ,5) 

CALL  TRISTF(1 ,2,5) 

CALL  TRI3TF( 2 ,3,5) 

CALL  TRI3TF(3,4,5) 

DO  91  1=1,6 
DO  91  J  =  1 ,15 

91  BS(I,J)=BS(I,J)/Afi 

DO  300  1=1 ,30 
DO  300  J= 1 ,30 
300  S1(I,J)=0.0 

AR1  =0.0 

CALL  NAXSTF(4 ,1,5) 

CALL  NAXSTF(1 ,2,5) 

CALL  NAX3TF(2 ,3,5) 

CALL  NAXSTF( 3 ,4,5) 

DO  310  1=1,6 
DO  310  J=1 ,30 

310  BS 1 ( I , J ) =  BS1 (I, J)/AR1 
DO  320  1=1,6 
DO  320  J=1 ,3 

320  BS 1 T ( I , J )  =  BS1 (I, J+12) 

DO  325  1=1,6 
DO  325  J=1 , 12 

325  BS1 (I, J+12)  =  BS1(I,J+15) 

DO  330  1=1,6 
DO  330  J  =  1 , 3 
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330  BS1 (1, J+24)  =  BS1T(IfJ) 

DO  340  1=1,3 

340  PIT  (1)  =  P1(I+12) 

DO  341  1=1,12 

341  PKl+12)  =  PI  (1+15) 

DO  342  1=1,3 

342  PI (1+24)  =  P IT (I ) 

DO  149  1=1,3 

DO  149  J=1 ,30 

149  S1T£M( 1 , J)  =  SI (1+12 , J ) 

DO  151  1=1,12 
DO  151  J=1 , 30 

151  SI (1+12 , J)  =  S1(1+15,J) 

DO  152  1=1,3 

DO  152  J=1 , 30 

152  SI (1+24 , J)  =  S1T£M(1,J) 

DO  153  1=1,3 

DO  153  J=1,30 

153  S1TErt(I, J)  =  S 1 ( J , 1+ 1 2 ) 

DO  154  J  =  1 , 12 

DO  154  1=1,30 

154  SI ( 1 , J+12 )  =  31(1, J+15) 

DO  155  1=1,3 

DO  155  J=1 , 30 

155  S1(  J , 1+24 )  =  S1TEM(I, J) 

DO  251  1=1,6 
DO  251  J=1 ,24 

251  TS(I, J)  =  0.0 

DO  252  1=1,3 
DO  252  J  =  1 ,4 
TS(I,I+( J-1 )*3)  =  0.250 

252  TS(I+3,I+12+(J-1)*3)  =  0.250 

DO  253  1=1,24 
DO  253  J= 1 , 24 
3 1 T ( I , J )  =  0.00 
DO  253  K=1 ,6 

253  S1T(1, J)  =  S 1 T ( 1 , J )  +  S1(I,24+K)*TS(K,J) 

DO  254  1=1,24 

DO  254  J=1 ,24 

254  S1(I,J)  =  S1(I,J)  +  S1T(I, J)  +  S 1 T ( J , I ) 

DO  255  1=1,24 

DO  255  J  = 1 , 6 
S2T(I , J)  =  0.0 
DO  255  K=1 ,6 

255  S2T(I,  J)  =  S2T(I , J)  +TS(K , I)*S1 (K+24 , J+24 ) 

DO  256  1=1,24 

DO  256  J=1 ,24 

S1T(I , J)  =  0.0 
DO  256  K= 1 , 6  '  ' 

256  S 1 T ( I ,  J )  =  S1T(I , J)  +  S2T(I,K)*TS(K, J) 
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DO  257  1=1,24 

DO  257  J  =1,24 

257  S1(I,J)  =S1(I,J)+S1T(1,J) 

RETURN 

END 

SUBROUTINE  SOLV 

COMMON /ELD AT A/ BETA ( 1000) ,EPR( 1 000) , PR (200) ,3H(200) , IX( 1 000 ,5) , 

IIP (200) , JP(200) , IS (200) , JS(200) , ALPHA ( 1000) , IT (200) , JT(200) , 

2ST (200 ) 

COMMON/BASiC/ACELZ , ANGVEL , ANGACC , THEF , VOL , NUMNP , NUMEL , NUMPC , NUMSC , 
1NUMST 

COMMON /N0NAX1 /SI (30,30) , PI (30) , THETA, BS1 (6,30) 

COMMON/NXDATA/NTP , NT  YPS , NTS , NT0T3 , GTS 1G ( 24 , 24 , 4 ) 

COMMON /SOLVE/ B( 162) , A( 1 62 , 81 ) , NUMBLK , MBAND 

MM=MBAND 

NN=8l 

NL=NN+1 

NH=NN+NN 

REWIND  1 

REWIND  2 

NB=0 

GO  TO  150 

Q*  *********************************** 

C  REDUCE  EQUATIONS  BY  BLOCKS 

(j*  *********************************** 

C 

C  1 .  SHIFT  BLOCK  OF  EQUATIONS 

C 

100  NB=NB+1 

DO  125  N=1 , NN 

NM=NN+N 

B(N)=B(NM) 

B(NM)=0 .00 
DO  125  M=1 ,MM 
A(N,M)=A(NM,M) 

125  A(NM,M)=0.00 

2.  READ  NEXT  BLOCK  OF  EQUATIONS  INTO  CORE 

IF ( NUMBLK . EQ . NB )  GO  TO  200 
150  READ (2)  (B(N) ,(A(N,M) ,M=1 ,MM) ,N=NL,NH) 

IF(NB.EQ.O)  GO  TO  100 

5.  REDUCE  BLOCK  OF  EQUATIONS 
C 

200  DO  300  N=1 ,NN 

IF(A(N, 1 ) .EQ.O.OO)  GO  TO  300 
B(N )=B(N)/A(N , 1 ) 

DO  275  L=2,MM 

IF(A(N,L) .EQ.O.OO)  GO  TO  275 
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C=A(N,L)/A(N , 1 ) 

I=N+L-1 

J=0 

DO  250  K=L ,MM 
J=J+1 

250  A( I , J)=A(I, J)-C*A(N,K) 

6(1)=B(I)-A(N,L)*B(N) 

A(N ,L)=C 
275  CONTINUE 
500  CONTINUE 

4.  WRITE  BLOCK  OF  REDUCED  EQUATIONS  ON  FORTRAN  UNIT  1 

IF(NUMBLK.EQ.NB)  GO  TO  400 

WRITE  (1)  (B(N),(A(N,M),M=2,MM),N=1,NN) 

GO  TO  100 

************************************ 

back-substitution 

ft*********************************** 

400  DO  450  M=1 ,NN 
N=NN+1-M 
DO  425  K=2,MM 
L=N+K-1 

425  B(N)=B(N)-A(N,K)*B(L) 

NM=N+NN 
B(NM) =B(N ) 

450  A(NM,NB)=B(N) 

NB=NB-1 

IF(NB.EQ.O)  GO  TO  500 
BACKSPACE  1 

READ  (1)  (B(N),(A(N,M),M=2,MM),N=1,NN) 

BACKSPACE  1 
GO  TO  400 

************************************ 

ORDER  FORMER  UNKNOWNS  IN  B  ARRAY 

************************************ 

500  K=0 

DO  600  NB= 1 ,NUMBLK 
DO  600  N=1,NN 
NM=N+NN 
K=K+1 

600  B(K)=A(NM,NB) 

************************************ 

WRITE  SOLUTION 

************************************ 

NN 12  =  3*NUMNP 
1500  FORMAT ( "  ",5110) 

WRITE(26)  (B(I) ,1=1 ,NN12) 

WRITE (26 )((IX(I,J) ,J=1 ,4) ,1=1 , NUMEL) 

MPRINT=0 


115 


DO  710  N= 1 , NUMNP 
IF(MPRINT.NE.O)  GO  TO  700 
WRITE  (6,2000) 

MPRINT=59 

700  MPRINT=MPRINT-1 

710  WRITE  (6,2001  )  W-,B(i»N-2)  ,B(J»M-1 )  ,B(3»N) 

2000  FORMAT  (IjHI  NODAL  POINT , 1 8X , 2HUR , 1 8X , 2HUZ , 1 8X , 2HUT ) 

2001  FORMAT  (I13.3E20.7) 

RETURN 

END 

SUBROUTINE  STIFF 
INTEGER  CODE 

COMMON/BASIC/ACELZ , ANG VEL , ANGACC , TREF , VOL , NUMNP , NUMEL , NUMPC , NUM3C , 
1NUMST 

C0MM0N/ELDATa/BETA(1000) ,EPR(1000) ,PR(200) ,Sri(200) ,IX(1 000,5) , 

IIP (200) ,  JP(200) , IS (200) , JS(200) , ALPHA ( 1000) , IT (200) , JT(200) , 
23T(200) 

COMMON /NPDATA/R( 1000) , CODE (1000) ,XR(1000) ,Z( 1000) ,XZ( 1000) , 

1NPNUM( 25 , 80 ) ,T( 1 000 ) ,XT( 1 000) 

COMMON /SOL VE/B( 1 62 ) ,A( 1 62 , 81 ) , NUMBLK,MBAND 
COMMON/NXDATA/NTP , NTTP3 , NTS , NTOTS , GTS 1G (24 , 24 , 4 ) 
COMMON/ANS4/FT(24,4) ,GT31U(24) ,GTS1UT(24,4) 
COMMON/ARG/RRR(5),ZZZ(5),RR(4),ZZ(4),S(15,15),P(15) ,TT(6), 

1H(6 , 15 ) ,CRZ(6 , 6) ,XI( 10 ) , ANGLE(4 ) ,SIG( 1 8) ,EPS( 18 ) ,N 
C0MM0N/N0NAXI/S1  (30,30)  ,P1  (30) , THETA, BS1  (6,30) 

COMMoN/PLANE/NPP 

C0MM0N/ANS2/GTP1 (24) ,G(24,24) ,GTS1 (24 , 24 ) ,GTS1GE(24 , 24 ) 
COMMON/DUM1/S1TEM(3,30) ,S1T(24,24),TS(6,24) 

DIMENSION  LM(4),S2(12,3),S3(3,12),34(3,3),S5(12,3),S6(12,12) 

C*  *********************************** 

C  INITIALIZATION 
REWIND  2 
REWIND  3 
NB=27 
ND=3*NB 
ND2=2»ND 
ST0P=0 . 

NUMBLK=0 
DO  100  N=1 ,ND2 
B(N )=0 . 00 
DO  100  M=1 ,ND 
100  A(N ,M)=0 . 00 
DO  50  1=1,24 
FT(I,NTP)  =  0.0 
GT31UT(I,NTP)  =  0. 0 
DO  50  J=1 ,24 

50  GTS1G(I, J,NTP)  =  0.0 

C*  *******»***************»***»»*»»»** 

C  FORM  STIFFNESS  MATRIX  IN  BLOCKS 

G«  *«»*«*««*»»»«»***««**«»************ 
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200  NUMBLK=NUMBLK+1 
NH=NB* (NUMBLK+1 ) 

NM=NH-NB 

NL=NM-NB+1 

KSHIFT=3*NL-3 

DO  340  N=1  ,NUMEL 

IF(IX(N,5) .LE.O)  GO  TO  340 

DO  210  1=1,4 

IF(1X(N,I) .LT.NL)  GO  TO  210 
IF(IX(N,I) .LE.NM)  GO  TO  220 
210  CONTINUE 
GO  TO  340 
220  CALL  QUAD 

IF(VOL.GT.O.)  GO  TO  230 
WRITE (6, 2000)  N 
ST0P= 1. 

230  IF( IX(N , 3) .EQ. IX(N , 4 ) )  GO  TO  300 
DO  231  H=1,3 

DO  231  JJ=1 , 3 

231  34(11, JJ)=S(II+12, JJ+12) 

CALL  SYMINV(S4,3) 

DO  232  11=1,12 
DO  232  JJ=1 ,3 

232  32(11, JJ)=S(Ii, JJ+12) 

DO  233  11=1,3 

DO  233  JJ=1 , 12 

233  33(H, JJ)=S(II+12,JJ) 

DO  240  1=1,12 

DO  240  J=1 , 3 
S5(I,J)=0.00 
DO  240  K=1 , 3 

240  35(1, J)  =  S5(I,J)  +  S2(I,K)  *  S4(K,J) 
DO  241  1=1,12 

DO  241  J = 1 , 12 
S6(I, J)=0.00 
DO  241  K=1 ,3 

241  S6(I, J)  =  36(1, J)  +  S5(I,K)  *  S3(K,J) 
DO  234  11=1,12 

DO  234  JJ=1 , 3 

234  P(1I)=P(II)-S5(II,JJ)*P( JJ+12) 

DO  235  11=1,12 

DO  235  JJ=1 , 12 

235  S(II , JJ)=S(II , JJ)-S6(II , JJ ) 

DO  259  1=1,24 

DO  259  J=1 , 24 
259  G(I, J)  =  0.0 
DO  260  K= 1 , 4 
DO  260  1=1,3 
G(K*3-3+I , 1*4-3 )  =1.0 
G(K*3-3+I, 1*4-2)  =  RRR(K) 
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G(K*3-3+I, 1*4-1)  =  ZZZ(K) 

260  G(K*3-3+I,I*4  )  =  ZZZ(K)  *RRR(K) 

DO  262  1=1 , 12 
DO  262  J  =  1 ,12 
26 2  G( 1+1 2 , J+1 2 )  =  G(I,J) 

NTP20  =  NT P+20 

WRITE (NTP20)  ((  CRZd.J)  ,J  =  1  ,6)  ,1=1  ,6) 

WR1TE(NTP20)((BS1(1,J) ,J=1 ,30) ,1=1 ,6) 

WR±TE(NTP20)((G(1,J) , J  =  1 , 24 ) , 1=1  , 24 ) 

DO  280  1=1,24 
GTP1 (I)=0.0 
DO  280  K=1 ,24 
GTSKI.K)  =  0.0 
GTP1 (i)=  GTP1(I)+  G(K , I) *P 1 (K ) 

DO  280  J=1 ,24 

280  GTSKI.K)  =  GTSKl.K)  +  G(J,I)  *  S1(J,K) 

WRITE(3)  ((GTS1 (I, J) ,J=1 ,24) ,1=1 ,24) 

DO  281  1=1 ,24 

FT(i , NTP)  =FT(I ,NTP)  +  GTP1(1) 

DO  281  J=1 ,24 
GTS 1GE( I , J )  =  0.0 
DO  281  K= 1 ,24 

281  GTS1GE(I, J)  =  GTS1GE(1 , J)+  GTSld.K)  *G(K,J) 

DO  282  1=1 ,24 
DO  282  J=1 ,24 

282  GTSIGd,  J.NTP)  =  GTS1G(I ,  J , NTP)  +  GTS1GE(I,J) 

C*  *********************************** 

C  ADD  ELEMENT  STIFFNESS  MATRIX  TO  BODY  STIFFNESS  MATRIX 

C*  *********************************** 
300  DO  310  1=1,4 

310  LM(I)=3*IX(N , I ) —3 
DO  330  1=1,4 
DO  330  K=1,3 
II=LM(I )+K-KSHIFT 
KK=3*l-3+K 
B(II)=B(II)+P(KK) 

DO  330  J=1 ,4 
DO  330  L=1 , 3 
J J=LM( J )+L-II+1 -KSrilFT 
LL=3*J-3+L 

IF(JJ.LE.O)  GO  TO  330 
TF(ND.GE.JJ)  GO  TO  320 
WRITE(6 ,2001 )  N 
STOP= 1 . 

GO  TO  3^0 

320  A(II,JJ)=A(II,JJ)+S(KK,LL) 

330  CONTINUE 
340  CONTINUE 

C*  *********************************** 

C  ADD  CONCENTRATED  FORCES 


118 


Cj  o 


c*  *********************************** 

DO  400  N=NL ,NM 
IF(N.GT.NUMNP)  GO  TO  500 
K=i*N-KSHIFT 
B(K)=B(K)+XT(N) 

B(K-1)=B(K-1)+XZ(N) 

400  B(K-2 )=B(K-2 )+XR(N ) 

************************************ 

ADD  PRESSURE  BOUNDARY  CONDITIONS 

C*  *********************************** 

500  IF(NUMPC.EQ.O)  GO  TO  600 
DO  540  L= 1 , NUMPC 
I=IP(L) 

J=JP(L) 

PP=PR(L)/6 . 

DR=(R( J)-R(I) )*PP 
DZ=(Z(I)-Z(J))*PP 
RX=2.*R(I)+R(J) 

ZX=R(I)+2.*R(J) 

11= 3*1 -KSH1FT-1 
JJ=i*J-KSHlFT-1 

IF  ( II .  LE .  0 .  OR .  II .  GT .  ND )  GO  TO  520 
SINA=0 . 

C03A=1 . 

510  B(II-1 )=B(II-1 ) +RX* ( COSA*DZ+SINA*DR ) 
GR=RX*(COSA*DZ+SINA*DR)*THETA/2 . 0 
FT ( 1 ,NTP)  =  FT( 1 , NTP)+GR 
FT(2 , NTP)  =  FT(2 , NTP)  +R(I)*GR 
FT(3,NTP)  =  FT(3,NTP)  +Z(I)*GR 
FT(4 ,NTP)  =  FT(4 , NTP)  +  R(I)*Z(I)*GR 
FT ( 1 4 , NTP)=  FT (14 ,NTP)  +R(I)*GR 
FT(13,NTP)  =  FT( 13, NTP)  +GR 
FT( 15 , NTP)  =  FT (15, NTP)+Z( I)*GR 
FT ( 1 6 , NTP )  =  FT ( 1 6 ,NTP)  +Z(I)*R(I)*GR 
B( II ) =B( II ) -RX* ( SINA*DZ-COSA*DR ) 

GZ=-RX* ( SINA*DZ-COSA*DR )  *THETA/2 . 0 
FT (5 , NTP)=FT(5 ,NTP)+GZ 
FT ( 6 , NTP)=FT(6 , NTP)+R(I)*GZ 
FT (7 ,NTP)=FT(7 ,NTP)+  Z(I)*GZ 
FT ( 8 , NT  P ) =FT ( 8 , NTP ) +Z ( I ) *R ( I ) *GZ 
FT (17, NTP)=FT( 1 7 , NTP)+GZ 
FT ( 1 8 ,NTP)=FT ( 1 8 ,NTP)+R (I )*GZ 
FT (19 ,NTP)=FT ( 1 9 ,NTP)+Z(I)*GZ 
FT (20 ,NTP)=FT(20 , NTP)+Z(I)*R(I)*GZ 

520  IF(JJ.LE.O.OR.JJ.GT.ND)  GO  TO  540 
3INA=0. 

COSA= 1 . 

530  B( JJ-1 )=B( JJ-1 )+ZX* ( COSA*DZ+SINA*DR) 

GR=  ZX  * ( C03A*DZ+SINA*DR )  *THETA/2.0 
FT ( 1 ,NTP)=FT(1 ,NTP)+GR 
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FT(2,NTP)=FT(2,NTP)+R(J)*GR 
FT(3 ,NTP)=FT(3 ,NTP)+Z( J )*GR 
FT(4,NTP)=FT(4,NTP)+Z( J)*R(J)*GR 
FT( 1 3 , NTP) =FT (13, NTP )+GR 
FT( 14,NTP)=FT( 14,NTP)+R(J)*GR 
FT ( 1 5 , NTP ) =FT (15,NTP)+Z(J) *GR 
FT ( 1 6 , NTP)=FT( 1 6 ,NTP)+Z( J )*R( J)*GR 
B ( J J ) =B( J J ) -ZX* ( SXNA*DZ-COSA*DR ) 

GZ=  -ZX* (SINA*DZ-C08A*DR )  *THETA/2.0 
FT(5 ,NTP)=FT(5 , NTP)+GZ 
FT(6 ,NTP)=FT(6 ,NTP)+R( J)*GZ 
FT (7 ,NTP)=FT(7 , NTP)  +Z(J)*GZ 
FT(8 ,NTP)=FT(8,NTP)+Z( J)*R( J)*GZ 
FT(17,NTP)=FT(17,NTP)+GZ 
FT( 1 8,NTP)=FT( 18 ,NTP)+R(J )*GZ 
FT (19 , NTP)=FT( 1 9 ,NTP)+Z( J ) *GZ 
FT(20,NTP)=FT(20,NTP)+Z( J)*R(J)*GZ 

540  CUNTINUE 

1100  FORMAT ( "  " , 12E10.3) 

C*  *********************************** 

0  ADD  SHEAR  BOUNDARY  CONDITIONS 

C*  *********************************** 

600  IF(NUMSC.EQ.O)  GO  TO  701 
DO  640  L=1,NUMSC 
I=IS(L) 

J=JS(L) 

SS=SH(L)/6 . 

DZ=(Z( I)-Z( J) )*SS 
DR=(R(J)-R(I))*SS 
RX=2.*R(I)+R(J) 

ZX=R(I)+2.*R(J) 

II=3*I-KSHIFT-1 

JJ=3*J-KSHIFT-1 

IF ( II . LE . 0 . OR . II . GT . ND )  GO  TO  620 
SINA=0 . 

COSA= 1 . 

610  B(II-1 )=B(lI-1 )+RX* ( SINA*DZ+C03A*DR ) 

GR=  RX*(SINA*DZ+COSA*DR)  *THETA/2.0 
FT ( 1 ,NTP)=FT(1 ,NTP)+GR 
FT(2, NTP) =FT (2, NTP)+R (I )*GR 
FT(3,NTP)  =FT(3 ,NTP)+Z(I)*GR 
FT (4 , NTP) =FT (4 ,NTP)+Z( I) *R(l )*GR 
FT(13,NTP)=FT(13,NTP)+GR 
FT ( 1 4 , NTP ) =FT ( 1 4 , NTP ) +R ( I ) *GR 
FT( 1 5 ,NTP)=FT ( 1 5 , NTP)+Z(I )*GR 
FT ( 1 6 ,NTP)=FT ( 1 6 , NTP)+Z( I)*R(I)*GR 
B( II )=B(II )-RX*(COSA*DZ-SINA*DR) 

GZ=  -RX* ( C0SA*DZ-3INA*DR )  *THETA/2.0 
FT(5,NTP)=FT(5,NTP)+GZ 
FT(6 ,NTP)=FT(6 ,NTP)+R(I)*GZ 
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FT (7 ,NTP)=FT(7 ,NTP)+Z(I)*GZ 
FT(8,NTP)=FT(8,NTP)+Z(I)*R(I)*GZ 
FT ( 1 7 , NTP ) =FT ( 1 7 , NTP ) +GZ 
FT(18,NTP)=FT(18,NTP)+R(I)*GZ 
FT( 1 9 ,NTP)=FT ( 1 9 ,NTP)+Z( i)*GZ 
FT (20 , NTP)=FT (20 , NTP)+Z( I ) *R ( I  )*GZ 

620  IF(JJ.LE.0.0R. JJ.GT.ND)  GO  TO  640 
SINA=0 . 

COSA=1 . 

630  B( JJ-1 )=B( JJ-1 )+ZX* ( 3INA*DZ+C0SA*DR ) 
GR=  ZX*(SINA*DZ+C0SA*DR)*THETA/2 . 0 
FT ( 1 , NTP)=FT( 1 , NTP)+GR 
FT ( 2 , NTP ) =FT ( 2 , NTP ) +R ( J ) *GR 
FT(3,NTP)=FT(3,NTP)+Z(J)*GR 
FT(4,NTP)=FT(4,NTP)+Z(J)*R(J)*GR 
FT ( 1 3  »NTP)=FT (13, NTP)+GR 
FT(14,NTP)=FT(14,NTP)+R(J)*GR 
FT ( 1 5 , NTP)=FT( 1 5 , NTP)+Z( J) *GR 
FT( 1 6 , NTP)=FT ( 1 6 ,NTP)+Z( J)*R( J)*GR 
B( J J ) =B( J J ) -ZX* ( C0BA*DZ-SINA*DR ) 

GZ=  -ZX* ( COSA*DZ-SINA*DR ) *THETA/2 . 0 
FT(5 , NTP)=FT (5 , NTP)+GZ 
FT(6,NTP)=FT(6,NTP)+R(J)*GZ 
FT (7 , NTP)=FT (7 ,NTP)+Z( J)*GZ 
FT(8,NTP)=FT(8,NTP)+Z(J)*R(J) *GZ 
FT (17 ,NTP)=FT( 1 7 ,NTP)+GZ 
FT(18,NTP)=FT(18,NTP)+R(J)*GZ 
FT( 1 9 ,NTP)=FT (19, NTP)+Z( J)*GZ 
FT(20,NTP)=FT(20,NTP)+Z(J)*R(J)*GZ 

640  CONTINUE 

701  IF(NUMST.EQ.O)  GO  TO  700 
DO  680  L= 1 , NUMST 
I=IT(L) 

J=JT(L) 

RT=ST(L)/6. 

RX=2.*R(I)+R(J) 

ZX=R(I )+2 . *R ( J) 

XX=SQRT ((R(J)-R(I))**2+(Z(J)-Z(I))**2) 

ii=3*i-kshift 

jj=3*j-kshift 

IF  ( II .  LE .  0 .  OR .  II .  GT .  ND )  GO  TO  670 
B(II)=B( II )+RT*RX*XX 
GT=RT*RX*XX*THETA/2 . 0 
FT (9 ,NTP)=FT(9 ,NTP)+GT 
FT( 10 , NTP)=FT( 10 , NTP)+R (I)*GT 
FT( 1 1 ,NTP)=FT( 1 1 , NTP)+Z(I)*GT 
FT( 1 2 ,NTP) =FT ( 1 2 , NTP)+Z( I)*R(I)*GT 
FT (21 ,NTP)=FT(21 ,NTP)+GT 
FT(22,NTP)=FT(22,NTP)+R(I)*GT 
FT(23  >  NTP )=FT (23 , NTP)+Z(I )*GT 
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FT(24,NTP)=FT(24,NTP)+Z(I )*R(I)*GT 
670  IF(JJ.LE.0.0K. JJ.GT.ND)  GO  TO  680 
6( JJ )=B( JJ )+RT*ZX*XX 
GT=RT*ZX*XX*THETA/2 .0 
FT(9,NTP)=FT(9>NTP)+GT 
FT(10,NTP)=FT(10,NTP)+R(J)*GT 
FT (11 , NTP)=FT ( 1 1 ,NTP)+Z( J)*GT 
FT( 12 , NIP )=FT ( 12 , NTP)+Z( J )*R( J )*GT 
FT (21 ,NTP)=FT(21 ,NTP)+GT 
FT(22 , NTP)=FT(22 , NTP)+R ( J )*GT 
FT(23,NTP)=FT(23,NTP)+Z(J)*GT 
FT(24 ,NTP)=FT(24,NTP)+Z( J)*R(J)*GT 
680  CONTINUE 

c  ADO  DISPLACEMENT  BOUNDARY  CONDITIONS 

700  DO  750  M=NL,NH 
IDM=0 

IF(M.GT .NUMNP)  GO  TO  750 
IF ( CODE (M) .GT .3)  GO  TO  751 
U=XR (M) 

N  =  3*M-2 -KSHiFT 
752  IF ( CODE (M) )  740,750,710 
710  IF ( CODE (M) .EQ. 1 )  GO  TO  720 
IF  (CODE(M) .EQ. 2)  GO  TO  740 
IF ( CODE (M) .EQ. j)  GO  TO  730 
GO  TO  740 

720  CALL  M0DIFY(ND2,N,U) 

CODE (M) = CODE (M)+IDM 
Gu  TO  750 

730  CALL  M0DIFY(ND2, N, U) 

740  U=XZ(M) 

N=N+1 


CALL  M0DIFY(ND2,N,U) 
CODE ( M) =  CODE ( M) +IDM 
GO  TO  750 
751  IDM=IDM+4 
U=XT(M) 

N=3*M-KSHIFT 


CALL  MODIFY (ND2,N,U) 
U=XR(M) 

N = 3  *M-2 -KSHIFT 
IF ( CODE (M) .EQ.4)  GO  TO  750 
CODE(M)=CODE(M)-4 
GO  TO  752 
750  CONTINUE 


C* 

C 

C* 


# 

* 


WRITE  BLOCK  OF  EQUATIONS  ON  FORTRAN  UNIT  AND  SHIFT  UP  LOWER  BLOCK 

*****#************j(j(J(j(j(j(j(j(j(#j(j(j(#j(j( 

WRITE  (2)  (B(N) ,(A(N,M) ,M=1 ,MBAND) ,N=1 ,ND) 
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DO  800  N=1 ,ND 

K=N+ND 

B(N)=B(K) 

B(K)=0.00 
DO  800  M= 1 , ND 
A(N,M)=A(K,M) 

800  A(K,M)=0 .00 

C*  a******#*#************############# 

C  CHECK  FOR  LAST  BLOCK 

C*  *********########****####»»«»»»»»»# 


IF(NM.LT.NUMNP)  GO  TO  200 
IF(STOP.NE.O. )  STOP 

2000  FORMAT  (27H  NEGATIVE  AREA  ELEMENT  NO.,l4) 

2001  FORMAT  (46H  BAND  WIDTH  EXCEEDS  ALLOWABLE  FOR  ELEMENT  NO.,14) 
RETURN 

END 

SUBROUTINE  STORE 
INTEGER  CODE 

COMMON /ANS4/  FT(24,4) ,GTS1U(24) ,GTS1UT(24,4) 

COMMON /NXDATA/NTP , NT IPS , NTS , NTOTS , GTS 1G ( 24 ,24,4) 

COMMON/NXMESH/ THETAN(4 ) ,NST(4) ,NUMS(4,5) ,NPC(8,8) 

COMMON/NPDATA/R( 1000) .CODE ( 1000 ) ,XR( 1000 ) ,Z( 1000 ) ,XZ( 1000 ) , 
1NPNUM(25 , 80 ) ,T( 1000 ) ,XT (1000 ) 

COMMON/SOLVE/B( 1o2) ,A(1o2,8l ) , NUMBLK.MBAND 
COMMON/GLBSEG/FI(24,8),FE(24,8),UC(24,8),SK(24,24,8) 

COMMON/BASIC/ ACELZ , ANGVEL , ANGACC , TREF , VOL , NUMNP , NUMEL , NUMPC , NUMSC , 
1NUMST 

COMMON/ANS2/LW(24) ,R1 (24,24) ,SK1 (24 , 24) , DUMM(24 , 24) 

DIMENSION  MW(24) 

JJ  =NST(NTP) 

DO  100  L=1 , JJ 
DO  50  1=1,24 
DO  50  J=1 ,24 
50  fl1(I,J)  =  0.0 
NS  =  NUMS(NTP.L) 

DO  110  KK  =  1,4 
NP1  =  NPC(NS ,KK) 

NP2  =  NPC(NS ,KK+4) 

DO  110  1=  1,3 

R1(3*(KK-1)+1  ,1*4-3  )  =  1.0 

R1 (3*(KK-1 )+I  ,1*4-2  )  =  R(NP1 ) 

R1 (3*(KK-1 )+I  ,1*4-1  )  =  Z(NP1 ) 

R1 (3*(KK-1 )+I  ,1*4  )  =  R(NP 1 )  *  Z(NP1 ) 

R1 (3* (KK— 1 )+I+12 , 1*4+9  )  =  1.0 
R1 (3*(KK-1 )+I+12 ,1*4+10)  =  R(NP2) 

R1(3*(KK-1)+I+12, 1*4+11)  =  Z(NP2) 

R1(3*(KK-1)+I+12, 1*4+12)  =  R(NP2)»  Z(NP2) 

UC(3*(KK-1 )+I ,NS)  =  B(3*NP1-3+I) 

UC(3*(KK-1 )+I+12,NS)  =  B(3*NP2-3+I) 

110  CONTINUE 
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CALL  MINV(R1 ,24,D1 ,LW,MW) 

WRITE(25)  ((H1(I, J),J=1 ,24) ,1=1 ,24) 

DO  115  1=1 ,24 
F£(I ,NS)=0 . 0 
FI(1,NS)=0.0 
DO  115  J=1 ,24 

FE(I,NS)  =  FE(I,NS)+  R1 ( J, I)*FT( J ,NTP) 

FI(I ,NS)  =  FI(I, NS)  +  R1 ( J , I)*GTS1UT ( J , NTP ) 

SKI (I , J)  =  0.00 
DO  115  K= 1 , 24 

SKI (I , J)  =  SKI (I, J)  +  R1 (K , I)  *  GTS1G(K , J,NTP) 

115  CONTINUE 

DO  120  1=1 ,24 

DO  120  J=1 ,24 

SK(I,J, NS)  =0.0 
DO  120  K= 1 , 24 

SK(I,J,NS)  =  SKU,J,NS)  +  SKI  (I  ,K)  *  R1(K,J) 

120  CONTINUE 

C  WRITE(6, 1000) (FE(I ,NS) ,1=1 ,24) 

WRITE(6, 1000) (FI (I, NS) ,1=1 ,24) 

WfiITE(6,1000)((SK(I,J,N3),J=1 ,24) ,1=1  ,24) 

WRITE(6,1000)(UCU,N3)  ,1=1  ,24) 

100  CONTINUE 

1000  FORMAT ( "  ",12E10.3) 

RETURN 

END 

SUBROUTINE  STRESS 
INTEGER  CODE 

COMMON/ANS4/FT(24,4) ,GTS1U(24)  ,GTS1UT(24,4) 

COMMON /BASIC/ ACELZ , ANGVEL , ANGACC , TREF , VOL , NUMNP , NUMEL , NUMPC , NUMSC , 
1NUMST 

COMMON /MATP/RO ( 6 ) ,E( 12, 16 ,6) ,EE( 16) ,AOFTS(6) 

COMM0N/NPDATA/R(1000) , CODEC  1 000 ) ,XR( 1000 ) ,Z( 1 000 ) ,XZ( 1 000) , 
1NPNUM(25 ,80) ,T(1000) , XT (1000) 

COMMON /ELDATA/BETAC 1000) ,EPR(1000) ,PR(200) ,SH(200) ,IX(1000,5) , 

IIP (200) ,JP(200) , IS (200) ,JS(200) ,ALPHA( 1000) , IT (200) ,JT(200) , 
2STC200) 

COMMON /ARG/RRR (5) ,ZZZ(5) ,RR(4) ,ZZ(4) ,S( 1 5 , 1 5 ) , P( 1 5 ) ,TT(6 ) , 

1H(6 , 15 ) ,CRZ(6 , 6) , XI (10) ,ANGLE(4) ,SIG(18) ,EPS(18),N 
COMMON /NONAXI/S 1(30,30) , PI (30) , THETA, BS1 (6 , 30) 

COMMON /NXDATA/NTP , NTYP3 , NTS , NTOTS , GTS  1 G ( 24 , 24 , 4 ) 

COMMON /NXMESH/THETAN(4) ,NST(4) ,NUMS(4,5) ,NPC(8,8) 

C0MM0N/ARG1/SIG1 (18) , EPS  1 ( 1 8) 

COMMON /SOLVE/B( 162) , A( 1 62 , 81 ) ,NUMBLK,MBAND 

COMMON/CONVRG/IDONE 

COMMON /PLANE/NPP 

C0MM0N/RESULT/BS(6 , 15) ,D(6 ,6) ,C( 6 , 6) , AR , BB( 6 , 9 ) , CNS ( 6 , 6 ) 

DIMENSION  LM(4) ,TP(6) ,TR(3,3) ,Q( 3 ) 

DIMENSION  QQ(3) 

COMMON/DUM1/S1TEM(3,30) ,GTS1(24,24) ,TS(6,24) 
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DIMENSION  P1 1(24) 

C*  •*»*•****************#*»*««****«*»• 

C  INITIALIZE 

C*  **»***********#**##*##############* 

REWIND  3 
XKE=0 . 

XPE=0 . 

MPRINT=0 
ERROR= . 005 
IDONE=1 

DO  200  N=1 , NUMEL 
IX (N , 5)=IABS(IX(N , 5) ) 

C  CALL  QUAD 

MTYPE=IAB3(IX(N , 5) ) 

DO  100  1=1 ,4 
11=3*1 
JJ=3*IX(N , I) 

P1(II-2)  =  B( JJ-2 ) 

PI(II-I)  =  B( JJ-1 ) 

P1(II  )  =  B(JJ  ) 

PI (11+10)  =  B( JJ-2) 

PI (11+11)  =  B( JJ-1 ) 

PKII+12)  =  B(JJ) 

P(II-2)=B(JJ-2) 

P(II-1 )=B( JJ-1 ) 

100  P(II)  =B( JJ) 

READ(3 ) ( (GTS1 (x, J) , J=1 ,24) ,1=1 ,24) 

DO  115  1=1 ,24 
GT31U(I)=0.0 
DO  115  J=1 ,24 

115  GTSIU(I)  =  GTS1U(I)+  GTS1 (I, J)*P1 (J) 

DO  116  1=1,24 

116  GT31UT(I,NTP)=GT31UT(I,NTP)  +  GTSIU(I) 

GO  TO  200 

DO  110  1=1,3 
110  Q(I)=P(I+12) 

DO  120  1=1,3 
DO  120  J=1 ,3 

120  Tft(I,J)=S(I+12,J+12) 

CALL  SYMINV(TR,3) 

DO  125  J=1 , 3 

QQ( J)=0 . 00 

DO  125  K= 1 ,12 

QQ( J)=QQ( J)+S( J+12 ,K)*P(K) 

125  CONTINUE 

DO  130  1=1,3 
P(I+12)=0.00 
DO  130  J=1 , 3 

130  P(I+12)=P(I+12)+TR(I,J)*(Q(J) -QQ ( J ) ) 

500  CONTINUE 
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c 

C  MATRIX  P  NOW  CONTAINS  15  DISPLACEMENTS  FOR  QUADRILATERAL  ELEMENT 
C 

C  CALCULATE  AVERAGE  STRAINS 
C 

DO  140  1=1,6 
EPS(±)=0.00 
DO  140  J=1 , 15 

140  EPS ( I ) =EPS ( I ) +BS ( I , J ) *P ( J ) 

C 

C  CALCULATE  AVERAGE  STRESSES 

C 

DO  151  1=1,6 
SlG(I)=0.00 
DO  151  J=1 ,6 

151  SIG(I)=SIG(I)+CRZ(I,J)*EPS(J) 

DO  152  1=1,6 
152  SIG(I)=SIG(I)-TT(I) 

C 

C  CALCULATE  STRAINS  IN  N-S-T  COORDINATES 

C 

DO  150  1=1  ,6 
EPS(I+6)=0.00 
DO  150  J=1  ,6 
DO  150  K=1 ,6 

150  EPS(I+6 )=EPS(I+6)+D(T , J)*C( J ,K)*EPS(K) 

C 

C  CALCULATE  STRESSES  IN  N-S-T  COORDIATES 

C 

DO  1 60  1=1,6 
SIG(I+6)=0.00 
DO  1 60  J=1 ,6 

160  SIG ( 1+6 ) =SIG ( 1+6 ) +CNS ( I , J ) *EPS ( J+6 ) 

DO  1 5 1  M=1 ,6 
P (M) =0 . 00 
DO  1 6 1  11=1,3 

IF( AOFTS(MTYPE) .EQ. 1 . )  P(M)=CNS(M,II)*EE(II+9) 

161  P(M)=P(M)+(T(N)-TREF)*CN3(M,II)*EE(II+9) 

DO  152  1=1,6 

162  SIG(I+6)=SIG(I+6)-P(I) 

C 

C 

DO  300  1=1 ,12 
300  EPS(I)=100.0*EPS(I) 

IF ( MPRINT . NE . 0 )  GO  TO  210 
WRITE(6,2000) 

WRITE(6 ,2002) 

MPRINT  =  19 

210  MPRINT=MPRINT-1 

WRITE (6 ,2001 )  N,RRR(5) ,ZZZ(5) ,(SIG(I) ,1=1,12) 
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WHlT£(6,200i)  T(N) , (EPS(I) ,1=1 , 12) 

200  CONTINUE 

2000  FORMAT ( 129H1  EL  R  Z  SIGMAR  SIGMAZ  3XGMAC  SIGMA 

1RZ  SiGMAZC  SIGMACR  SIGMAN  SIGMAS  SIGMAT  SIGMANS  SiGMAST 

2  SIGMA IN ) 

2001  F0RMAT(1H0,I5, 1X.2F7-4, 12F9.0) 

2002  FORMAT ( 1 28H0  TEMPERATURE  EPSR  EPSZ  EP3C  EP3R 

1Z  EPSZC  EPSCR  EP3N  EPSS  EPST  EP3NS  EPSST 

2  EPSTN) 

2003  F0RMAT(6X,F13.0,2X, 12F9.5) 

RETURN 

END 

SUBROUTINE  SYMINV(A, NMAX) 

DIMENSION  A(NMAX,NMAX) 

DO  300  N=1 ,NMAX 
D=A(N ,N) 

DO  100  J=1 ,NMAX 
100  A(N, J)=-A(N, J)/D 
DO  210  1=1 , NMAX 
IF(N.EQ.I)  GO  TO  210 
DO  200  J= 1 , NMAX 

IF(N.NE.J)  A(I,J)=A(I,J)+A(I,N)*A(N,J) 

200  CONTINUE 
210  A(I,N)=A(I,N)/D 
300  A(N,N)=1 .00/D 
RETURN 
END 

SUBROUTINE  TEMP(R,Z,T) 

COMMON/SOLVE/X(4428) ,Y(4428) ,TEM(4428) ,NUMTC ,MBAND 
DIMENSION  SMALL (20) , ISM (20) 

C*  *********************************** 

INITIALIZE 

XXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXXXXX 

J=1 

JMaX=16 

IF (NUMTC . LT . JMAX)  JMAX=NUMTC 
DO  10  1=1, JMAX 
SMALL (I)=0. 

10  ISM( I )=0 

QXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C  FIND  THE  JMAX  CLOSEST  POINTS 

QX  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

DO  50  1=1, NUMTC 
DSQ=(X(1)-R)**2+(Y(I)-Z)**2 
IF(DSQ.GT . . IE-4)  GO  TO  20 
T=T£M(I) 

RETURN 

20  IF(I.EQ.I)  SMALL ( 1 )=DSQ 
IFQ.EQ.1)  ISM(  1  )=  1 
IF(I.EQ.I)  GO  TO  50 
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IF ( SMALL ( J ) .  LE . DSQ . AND . J . LT . JMAX )  SMALL(J+1 )=D8Q 
IF ( SMALL ( J ) . LE . DSQ . AND . J . LT . JMAX )  ISM(J+1 )=l 
IF ( SMALL (J ) .LE . DSQ)  GO  TO  40 
DO  30  K=1 , J 
JB=J-K  +1 

IF( JB.EQ.O)  GO  TO  40 
SMALL(JB+1 )=SMALL(JB) 

ISM( JB+1 )=ISM( JB) 

SMALL (JB)=DSQ 
ISM( JB) =1 

IF(JB.EQ.I)  GO  TO  40 
IF ( SMALL (JB-1 ) .LE.DSQ)  GO  TO  40 
30  CONTINUE 

40  IF(J.LT.JMAX)  J=J+1 
50  CONTINUE 

C*  *********************************** 

C  FIND  THE  THIRD  TEMPERATURE  POINT  BY  AREA  TEST 

C*  *********************************** 

JCHK=JMAX-2 

J=0 

I 1 =ISM( 1 ) 

I2=ISM(2) 

60  I3=ISM( J+3) 

AREA=.50*(Y(I1)*X(I3)-Y(I3)*X(n)+Y(I3)*X(I2)-Y(I2)*X(I3)+ 

1  Y(I2)*X(I1)-Y(I1)*X(I2)) 

D1=(X(12)-X(I1))**2+(Y(I2)-Y(I1))**2 

IF  D1  IS  APPROXIMATELY  0.  IT  IS  ASSUMED  THAT  THERE  EXISTS  A 
DUPLICATION  OF  INPUT 
IF(D1 .GT..1E-3)  GO  TO  70 
12=13 
J=J+1 
GO  TO  60 

70  IF(AREA**2.GT. . 1 *D1 “SMALL (1 ))  GO  TO  80 
J=J+1 

IF( J.LT.JCHK)  GO  TO  60 
WRITE(6 , 2000 )  II, 12, 13, J 
T=TEM(I1 ) 

RETURN 

C*  *********************************** 

C  FIND  TEMPERATURE  INTERCEPT 

C*  *********************************** 

80  DETA=Y(I1)*(TEM(I3)-TEM(I2))+Y(I2)*(TEM(I1)-TEM(I3)) 

1  +Y(I3)*(TEM(I2)-TEM(I1)) 

DETB=X(I1)*(TEM(I2)-TEM(I3))+X(I2)*(T£M(I3)-TEM(I1)) 

1  +X(I3)*(TEM(I1)-TEM(I2)) 

DETC=TEM(I1 )* (X (12 ) *Y( 13) -X ( 13) *Y ( 12 ))+TEM( 12)* (X( 13 )*Y( II ) -X ( I 1 )* 
1Y(I3))+TEM(I3)*(X(II)*Y(I2)-X(I2)*Y(II)) 

T  =  (DETA*R+DETB* Z+DETC ) / ( 2 . *  AREA ) 

2000  FORMAT  (28H  ERROR  IN  TEMPERATURE  INPUT, 5H  I1=I4,5H  12=14, 

15H  13=14, 4H  J=I4) 
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RETURN 

END 

SUBROUTINE  TEM2(NUMNP) 

INTEGER  CODE 

COMMON/NPDATA/R( 1000) ,CODE( 1 000) ,XR( 1000) ,Z( 1000) ,XZ(1 000) , 
1NPNUM(25,80) ,T( 1000) , XT (1000) 

READ(5 , 1000 )  T CON ST 
DO  100  N  =  1 , NUMNP 
100  T(N )=TCONST 
1000  FORMAT (FI  0.0) 

RETURN 

END 

SUBROUTINE  TRlSTF  (II.JJ.KK) 

INTEGER  CODE 

C0MM0N/MATP/R0(6) ,E(12, 1 6 , 6 ) , EE ( 1 6 ) ,A0FTS(6) 

COMMON /BASIC/ ACELZ , ANGVEL , ANGACC , TREF , VOL , NUMNP , NUMEL , NUMPC , NUMSC , 
1NUMST 

COMMON/ ARG/ RRR(5 ) ,ZZZ(5 ) ,RR(4),ZZ(4),S(15,15),P(15),TT(6), 

1 H ( 6 , 1 5 ) , CRZ ( 6 , 6 ) , XI ( 1 0 ) , ANGLE ( 4 ) , SIG (1 8 ) , EPS (1 8 ) ,  N 
COMMON/NPDATA/ROOOO) ,CODE( 1000), XR( 1000) ,Z( 1000) ,XZ( 1000) , 
1NPNUM(25,80) ,T(1 000) ,XT (1000) 

COMMON /ELD AT A/ BET A ( 1 000 ) ,  EPR( 1 000 ) , PR (200 ) ,SH(200 ) , IX( 1 000 ,5) , 

1 1P (200 ) , JP (200 ) , IS (200 ) ,  JS(200 ) , ALPHA ( 1000 ) , IT (200 ) , JT(200 ) , 
2ST(200 ) 

COMMON /NONAXI/S 1 ( 30 , 30 ) , PI  ( 30 )  , THETA , BS 1 ( 6 , 30 ) 
COMMON/RESULT/BS(6,15),D(6,6),C(6,6),AR,BB(6,9),CNS(6,6) 

DIMENSION  B1A(6,9),B1B(6,9),B2A(6,9),B2B(6,9),B3A(6,9),B3B(6,9) 
DIMENSION  B1(6,9),B2(6,9),B3(6,9),F(6,9),G(9,6),V(9,9) 

DIMENSION  BF(3) ,BFR(3) ,BFZ(3) , TP(9) ,B(9 ,9 ) 

MTYPE=IABS(IX(N , 5) ) 

RR(1 )=RRR(Ii) 

RR(2)=RRR(JJ) 

RR(3)=RRR(KK) 

ZZ(1 )=ZZZ(II) 

ZZ(2)=ZZZ( JJ ) 

ZZ(3)=ZZZ(KK) 

CALL  INTER 
VOL=VOL+XI( 1  ) 

C0MM=RR(2)*(ZZ(j)-ZZ(1 ))+RR(1 )*(ZZ(2)-ZZ(3))+RR(3)*(ZZ(1 )-ZZ(2)) 

DO  10  1=1,6 
DO  10  J=1 ,9 
B1 (I , J)=0 . 00 
B2(I, J)=0.00 
10  B3(I,J)=0.00 

C  FILL  B1  MATRIX-CONSTANT  TERMS 
B1(1  ,1)  =  (ZZ(2)-ZZ(3))/C0MM 
B1(1 ,4)=(ZZ(3)-ZZ(1 ))/COMM 
B1(1,7)=(ZZ(1)-ZZ(2))/C0MM 
B1 (i , 1 )=B1 (1,1) 

B1 (3,4)=B1 (1,4) 
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B1 ( j ,7)=B1 (1,7) 

B1  (2,2)=(RR(3) -RR ( 2 ) ) / COMM 
B1(2,5)=(RR(1)-RR(i))/COMM 
B1  (2,8)=(RR(2)-RR(1 ))/COMM 
B1 (4 , 1 )=B1 (2,2) 

B1 (4 , 4)=B1 (2 , 5) 

B1  (4 , 7)=B1 (2,8) 

81(4 , 2)=B1 (1,1) 

B 1  ( 4 , 5 ) =B 1(1 ,4) 

B1 (4 ,8)=B1 (1,7) 

B1 (5 , 3)=B1 (4,1) 

B1  (5 , 6 ) =B 1 (4,4) 

B1 (5 ,9)=B1 (4,7) 

C  FILL  B2  MATRIX-1 /R  TERMS 

B2(3, 1 )=(1/COMM)*((ZZ(3)-ZZ(2))*RR(2)+(RR(2)-RR(3))*ZZ(2)) 
B2(3,4)=(1/COMM)*((ZZ(1)-ZZ(3))*RR(3)-(RR(1)-RR(3))*ZZ(5)) 
B2(3 , 7)= ( 1/COMM)*( (ZZ(2)-ZZ( 1 ) )*RR( 1 )+(RR( 1 )-RR(2) )*ZZ( 1 ) ) 
B2(o,3)=-B2(3,1) 

B2(6 , 6)=-B2(3 , 4) 

B2(6 , 9)=-B2(3 , 7 ) 

C  FILL  B3  MATRIX-Z/R  TERMS 

B3(3,1)=(RR(3)-RR(2))/COMM 
B3(3,4)=(RR(1)-RR(3))/COMM 
B3(3,7)=(RR(2)-RR(1 ))/COMM 
B3(6,3)=(RR(2) -RR ( 3 ) ) / COMM 
B3(6,6)=(RR(3)-RR(1))/COi'1M 
B3(6,9)=(RR(1 )-RR(2))/COMM 
AR=AR+XI( 1 ) 

DO  80  1=1,6 
DO  80  J=1 ,9 

80  BB(I, J)=B1 (I,J)*XI(1 )+B2(I,J)*XI(2)+B3(I,J)*XI(4) 

DO  81  K=1 ,6 

DO  81  1=1,3 

BS (K , 3*J J-3+I )=BB(K, 1+3 ) +BS (K , 3*J J-3+I ) 
B3(K,3*II-3+I)=BB(K,I)+BS(K,  3*H-3+I) 

81  BS(K , 3*KK-3+I ) =BB(K , 1+6 )+BS(K , 3*KK-3+I ) 

DO  220  1=1,6 

DO  220  J=1 ,9 

B1A(I, J)=B1 (I, J)*XI(1 )+B2(I,J)*XI(2)+B3(I, J)*XI(4) 
B2A(I,J)=B1(I,J)*XI(2)+B2(I,J)*XI(3)+B3(I, J)*XI(5) 
B3A(I,J)=B1(I,J)*XI(4)+B2(I,J)*XI(5)+B3(I,J)*XI(6) 

220  CONTINUE 
DO  200  1=1,6 
DO  200  K=1 ,9 
B1B(I,K)=0 .0 
B2B(I,K)=0.0 
B3B(I,K)=0.0 
DO  200  J=1 ,6 

B1B( I ,K)=B1B(I,K)+CRZ(I , J)*B1A( J ,K) 
B2B(I,K)=B2B(I,K)+CRZ(I, J)*B2A(J,K) 
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o  c: 


BiB(I,K)=B3B(I,K)+CRZ(I,J)*BiA(J,K) 

200  CONTINUE 
DO  230  1=1,9 
DO  230  K=1 ,9 
B(I,K)=0.0 
Do  230  J=1 ,6 

B(I,K)=B(I,K)+B1(J,I)*B1B(J,K)+B2(J,I)*B2B(J,K)+B3(J,1)*B38(J,K) 
230  CONTINUE 

ASSEMBLE  QUADRILATERAL  STIFFNESS  MATRIX,  S,  FROM  TRIANGULAR 
STIFFNESS  MATRIX,  B. 

IIM=3*ll-3 
JJM=3*JJ-3 
KKM=3*KK-3 
DO  120  1=1,3 
DO  120  J=1 ,3 

S ( IIM+I , IIM+J ) =B( I  ,J  )+S(IIM+I , IiM+J ) 

3 ( IIM+I , JJM+J ) =B( I  , J+3 )+3 ( IIM+I , JJM+J ) 

3 ( IlM+I , KKM+J )=B( I  , J+6 ) +S ( IIM+I , KKM+J ) 

S ( J JM+I , IIM+J ) =B ( 1+3 , J  ) +S ( J JM+I , IIM+J ) 

3 ( J JM+I , JJM+J ) =B( 1+3 , J+3 )+3 ( J JM+I , JJM+J ) 

S ( J JM+I , KKM+J )  =  B ( 1+3 , J+6 ) +S ( J JM+I , KKM+J ) 

3  ( KKM+i ,  UM+J )  =B ( 1+6 ,  J  )  +S  ( KKM+I ,  IiM+J ) 

S (KKM+I , JJM+J ) =B( 1+6 , J+3 )+S (KKM+I , JJM+J ) 

3 ( KKM+I , KKM+J ) =3 ( 1+6 , J+6 ) +S (KKM+I , KKM+J ) 

120  CONTINUE 

C  ASSEMBLE  BODY  FORCES  MATRIX 

BF( 1 )= ( ZZ( 3)*RR(2)-RR(3 )*ZZ(2) ) /COMM 
BF(2)= (ZZ( 1 )*RR(3)-RR( 1 )*ZZ( 3) ) /COMM 
BF(3 )= (ZZ(2 )*RR( 1 )-RR(2 )*ZZ( 1 ) ) /COMM 
BFR(1 ) = ( ZZ ( 2 ) -ZZ ( 3 ) ) / COMM 
BFR(2 )= ( ZZ( 3 )-ZZ( 1 ) )/COMM 
BFR(3)= ( ZZ ( 1 )-ZZ(2) )/COMM 
BFZ( 1 )= (RR(3 )-RR(2) )/COMM 
BFZ(2)= (RR( 1 )-RR(3) )/COMM 
BFZ(3)=(RR(2)-RR(1 ))/COMM 
C  BODY  FORCE  IN  Z-DIRECTION 

COMM=-ACELZ*RO ( MTYPE ) 

DO  140  1=1,3 
IIK=3*I-1 

140  TP(1IK)=C0MM*( BF(I)*XI(  1  )+BFR(I )*XI(7 )+BFZ( I)*XI(8) ) 

C  BODY  FORCE  IN  R-DIRECTION 

C0MM=ANGVEL**2*R0 (MTYPE ) 

DO  150  1=1,3 
L=3*I-2 

150  TP(L)=COMM*( BF(I )*XI(7 )+BFR(I)*XI(9 )+BFZ(I)*XI( 10)) 

C  BODY  FORCES  IN  TANG.  DIRECTION 

COMM=-ANGACC*RO (MTYPE) 

DO  160  1=1,3 
IIM=3*I 

160  TP ( IIM) =COMM* ( BF ( I ) *XI ( 7 )+BFR ( I ) *XI ( 9 )+BFZ ( I ) *XI (10)) 
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C  ADD  THERMAL  EFFECTS 
DO  1 61  J=1 ,9 
DO  1 6 1  K=1 ,6 

161  TP(J)=TP(J)+(XI(1)*B1(K,J)+XI(2)*B2(K,J) 

1+Xl(4)*B3(K,J))*TT(K) 

C  REARRANGE  TP  INTO  P -MATRIX, THE  BODY  FORCES  MATRIX 
K=3*iI-2 
L=3*JJ-2 
M=3*KK-2 
DO  170  1=1,3 
J=I-1 

P1(K+J  )  =  P1(K+J  )+TP( I  )*THETA/2 . 0  +TP( 1+6 )*THETA/4 . 
PUK+J+15)  =  PI  (K+J+1 5 )+TP( I  )*THETA/2.0  +TP( 1+6 )*THETA/4 . 
P1(L+J  )  =  P1(L+J  )+TP( 1+3) *THETA/2 . 0  +TP( 1+6 )*TH£TA/4 . 
P1(L+J+15)  =  PI (L+J+1 5 )+TP(I+3)*THETA/2 . 0  +TP( 1+6 )*THETA/4 . 
P 1 (M+J  )  =  P 1 (M+J  )+TP( 1+6 )*TH£Ta/2 . 0 
PKM+J+15)  =  PI  (M+J+15)+TP(I+6)*THETA/2.0 
P(K+J)=P(K+J)+TP(I) 

P(L+J)=P(L+J)+TP(I+3) 

170  P(M+J)=P(M+J)+TP(i+6) 

RETURN 

END 

SUBROUTINE  XMODFY(U.N) 

C0MM0N/NX30LV/SK  (132,24),R1  (132),NSZF 
N6AND=24 
DO  10  M=2 , NBAND 
K=N-M+1 

IF(K.LE.O)  GO  TO  5 
R1  (K)=R1 (K)-SK(K,M)*U 
SK(K,M)=0 . 

5  K=N+M-1 

IF(NSZF.LT.K)  GO  TO  10 
R1(K)=R1(K) -SK(N ,M) * J 
SK(N ,M)=0 . 

10  CONTINUE 
SK(N , 1 )=  1 . 

R1(N)=U 

RETURN 

END 

SUBROUTINE  XSOLVE 

COMMON /NX30LV/3K  ( 132,24) ,R1  (132),NSZF 
NBAND=24 
DO  300  N= 1 , N3ZF 
I=N 

DO  290  L=2, NBAND 
1=1+1 

IF(SK(N ,L) )  240,290,240 
240  AC=SK(N,L)/SK(N,1) 

J=0 

DO  270  K=L, NBAND 


0 

0 

0 

0 
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J=J+1 

IF(SK(N,K))  260,270,260 
260  SK(I,J)=SK(I,J)-AC«SK(N,K) 
270  CONTINUE 
280  SK(N,L)=AC 
C 

R1(I)=R1(I)-AC»R1(N) 

290  CONTINUE 

800  R1(N)=R1(N)/SK(N,1) 

C 

N=NSZF 
850  N=N-1 

IF(N )  500,500,360 
360  L=N 

DO  400  K=2 , NBAND 
L=L+1 

IF(SK(N ,K) )  370,400,370 
870  R1 (N)=R1 (N)-SK(N ,K)*R1 (L) 
400  CONTINUE 
GO  TO  350 
C 

500  RETURN 
END 
END 
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USER  EVALUATION  OF  REPORT 


Please  take  a  few  minutes  to  answer  the  questions  below;  tear  out 
this  sheet  and  return  it  to  Director,  US  Army  Ballistic  Research 
Laboratory,  ARRADCOM,  ATTN:  DRDAR-TSB,  Aberdeen  Proving  Ground, 
Maryland  21005.  Your  comments  will  provide  us  with  information 
for  improving  future  reports. 

1 .  BRL  Report  Number 

2.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related 
project,  or  other  area  of  interest  for  which  report  will  be  used.) 


3.  How,  specifically,  is  the  report  being  used?  (Information 
source,  design  data  or  procedure,  management  procedure,  source  of 
ideas,  etc.) _ 


4.  Has  the  information  in  this  report  led  to  any  quantitative 
savings  as  far  as  man-hours/contract  dollars  saved,  operating  costs 
avoided,  efficiencies  achieved,  etc.?  If  so,  please  elaborate. 


5.  General  Comments  (Indicate  what  you  think  should  be  changed  to 
make  this  report  and  future  reports  of  this  type  more  responsive 
to  your  needs,  more  usable,  improve  readability,  etc.) 


6.  If  you  would  like  to  be  contacted  by  the  personnel  who  prepared 
this  report  to  raise  specific  questions  or  discuss  the  topic, 
please  fill  in  the  following  information. 
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Telephone  Number: 
Organization  Address: 


